Assessment of reactivities with explicit and implicit solvent models: QM/MM and gas-phase evaluation of three different Ag-catalysed furan ring formation routes

Abstract

QM/MM molecular dynamics and static DFT calculations have been performed to evaluate free energy barriers and reaction free energies of a representative silver-catalyzed reaction in DMF. The mechanism developed in a previous work revealed a favorable intramolecular C–O coupling between terminal alkyne and β-ketoester moieties to yield a furan ring and three possible pathways were scrutinized within the framework of the SMD implicit solvation model. In this study we set out to compare the effect of implicit and explicit solvation on the three possible variations of the furan ring-closure step. The three pathways feature different charge states with bonding topologies characteristic of many silver- and copper-catalyzed reaction steps; hence they can serve as a blueprint for assessing the effects of solvents in a wider set of reactions. Comparison of the results showed that both methodologies could unequivocally determine the most favorable as well as the least likely pathway. Further analysis of the trajectories obtained from the QM/MM simulations indicated neither direct solvent participation in the reaction nor any site-specific interaction of the solvent with the reactant despite the fact that pairwise interactions between the solutes and the highly polar solvent molecules are significant. These insights point to a sufficiently mobile, fluctuating solvent shell which can be efficiently substituted by implicit solvent models at a fraction of computational costs.

---

Introduction

Collaboration between experiment and theory to understand and improve catalytic processes has become very efficient and productive over the last few decades.^1,2 An important requirement to achieve accurate theoretical rationalizations and predictions for experiment is to employ an adequate simulation model, i.e. one involving those atoms and interactions which play a decisive role in the mechanism. This requirement implies sufficiently accurate methods to calculate forces and sufficiently detailed molecular models. A particularly important issue is the solvent model employed in the calculations.^3,4 Implicit solvent models treat solvents as a continuous medium surrounding the solute, whereas explicit models take into account the movements and effects of the actual solvent molecules within a given region around the solute molecules.⁵ While implicit models speed up calculations enormously, many potentially important effects are neglected or averaged out.⁶ Clearly, when solvent molecules play a role not only in solvation but also in the reaction, explicit solvent models are the proper choice.⁷ However, when stoichiometry suggests no reactive solvent participation the choice is not straightforward: solvent molecules may have an effect on critical steps (steps determining the rate or selectivity) which cannot be captured by implicit models.⁸ In a recent study we have theoretically studied the Ag-catalysed oxidative furan formation from terminal alkyne precursors and β-ketoesters⁹ using the implicit solvent model of the nonreactive solvent DMF.¹⁰ The field of silver-catalysed coupling reactions has evolved notably recently and this reaction is a very good example of the high efficiency and versatility of Ag catalysis.¹¹ A crucial step of this synthetic route is the furan-ring formation. The ring-closure requires large atomic displacements which can influence the solvation shell. In the present study this step is evaluated using implicit and explicit solvent models to compare the performance of the two approaches. Scheme 1 displays the reaction routes probed by the calculations. An important aspect of these routes which makes them ideal for the present study is that they have different charge states: positive, neutral and negative, hence they represent the typical solute–solvent situations occurring in practice. As similar, analogous situations occur in many other silver- and copper-catalysed reactions,¹² our model systems can serve as representative examples for a larger number of reactions. Indeed, as Scheme 2 shows the selected analogous reactions also feature these metal ions bonded to a carbon atom and they are not shielded by their ligands, i.e. they are quite exposed to the solvent molecules. Such a pattern is important to consider when modelling the solvation effects. In addition, the oxidation state of the metal ions is the same in all cases. DMF is a prototype of polar, aprotic solvents (such as DMSO, acetonitrile, acetone, etc.) and can serve as a model for this family of solvents.¹³ As shown in this study, the essential result is that for this type of reactions the implicit solvation model is sufficient to rationalize the mechanism in terms of activation barriers and reaction free energies.

Scheme 1 Non-catalysed and silver-catalysed ring closure pathways studied in this work.

Scheme 2 Selected noble-metal catalyzed coupling reactions and a crucial transition state structure from the corresponding computed mechanistic analysis. (1) Ref. 12a; (2) ref. 12b and c; (3) ref. 12c and e; (4) ref. 12d; (5) ref. 12e.

Computational details

QM/MM calculations have been performed to explore the free energy profile of the reactions explicitly taking into account the effect of the solvent DMF molecules. The CP2K package has been used for these calculations.¹⁴ Reactants 1r, 2r and 3r (Scheme 1) were inserted into a periodic box with dimensions of 24.403 Å × 24.403 Å × 24.403 Å which hosted also 112 DMF molecules. The densities are 0.962 g cm⁻³ (1r) and 0.974 g cm⁻³ (2r, 3r). For each reaction the solute molecule has been described by DFT employing the PBE functional¹⁵ with the D3 dispersion corrections¹⁶ (PBE+D3). For the valence orbitals the short-range molecularly optimized double-ζ basis sets augmented with polarization functions were used.¹⁷ The electronic charge densities have been expanded in a plane-wave basis set with a cut-off of 400 Ry. The GTH pseudopotentials optimized to the selected functionals have been applied to describe the interactions between the ionic core and the valence electrons.¹⁸ The QM description of the solute has been restricted to a smaller cubic box with an edge of 20 Å. The flexible solvent molecules were treated with the CHARMM general force field¹⁹ and non-bonded parameters for all atoms have been taken from this force field. The electrostatic coupling between the QM and MM parts has been accounted for by the method of Laino et al.²⁰ The molecular dynamics simulations have been performed under NVT conditions using the velocity rescaling thermostat of Bussi et al.²¹ The time step was 0.5 fs. Blue moon sampling with thermodynamic integration (TI) has been done to obtain the free energy profiles of the reactions.²² The following protocol has been followed: after a 25 ps equilibration with the classical force field, where the solute molecule was fixed, a 10 ps QM/MM equilibration was done. Then an exploratory slow-growth simulation was performed with the distance between the reacting C and O atoms as the reaction coordinate. In these simulations the spring constant was 0.02 Ha Bohr⁻² for reactions 1 and 2, whereas 0.025 Ha Bohr⁻² for reaction 3. The target speed for decreasing the Ag–O distance was −0.0005 Å fs⁻¹ for reactions 1 and 2, and −0.0003 Å fs⁻¹ for reaction 3. The reaction coordinate spans the interval of 4.5–1.4 Å. Configurations from this trajectory have been selected for the TI-s where the C–O distance was fixed stepwise. The strategy for selection was the following: first configurations of regular intervals (250 fs for reactions 1 and 2, 150 fs for reactions 3) have been taken, then if it was necessary, additional configurations with suitable Ag–O distances have been selected for the TI. In this way 16, 13 and 15 different reaction coordinate points have been sampled for reactions 1, 2 and 3, respectively. Following a 2 ps equilibration the production runs were carried out for at least 5 ps. The uncertainties of the calculated free energies have also been estimated and it is found that they are always less than 0.7 kcal mol⁻¹.²³ The simulations have been performed for the experimental temperature of 353 K. For the visualization of the molecular dynamics results the TRAVIS program was used.²⁴

Calculations on the solute molecules have also been performed where the solvent effects were included with the SMD implicit solvent model.²⁵ These calculations have been done with the Gaussian 09 program²⁶ using both PBE+D3 and M06 functionals.²⁷ The M06 functional is considered to be more accurate for energetics,^28,29 and earlier calculations have employed it successfully.¹⁰ The 6-31G* basis set has been selected which corresponds in size and flexibility to the basis set employed in the QM/MM calculations. For Ag the valence electrons have been described by the LANL2DZ basis set augmented with a set of f-type polarization functions taken from the cc-pVDZ-PP basis set. The LANL2 effective core potential has been employed to substitute the ionic core of Ag. The solvent corrections have been included into the optimization and frequency calculations self-consistently. The free energy values are calculated assuming the validity of the ideal gas–harmonic oscillator–rigid rotator model as implemented in the Gaussian package. Note that this model to estimate entropy contribution (especially translational entropy) assumes that the degrees of freedom of the solute in gas phase and solution phase are very similar. This assumption is often proved to be inaccurate although in many other cases it works very well.³ The transition states (TS-s) have been determined by the TS optimization routine of the Gaussian package with the solvent corrections included into the optimization. All the optimization calculations have been verified by the frequency calculations (no imaginary frequency for the stable intermediates and a single imaginary frequency for the TS-s). Combined IRC and optimization calculations have been carried out to verify that the TS-s connect the corresponding intermediate states.

Results and discussion

The reaction steps featured in Scheme 1 are characterized energetically by the activation free energies and the reaction free energies collected in Tables 1 and 2, respectively. The three different mechanisms feature different charge states: overall negative, neutral and positive for reactions 1, 2, and 3, respectively. We note that these steps have been evaluated as possible reaction steps leading to the furan ring formation in ref. 10 already. It has also been shown that the most favourable pathway features reaction 2, whereas the other reactions (1 and 3) could be excluded. In particular reaction 1, where there is no silver participation, is highly unfavourable stressing the importance of the catalytic effect of Ag in this reaction. The principal aim of this study is to investigate the role of solvation in the reaction mechanism. To this end we performed QM/MM molecular dynamics simulations to explicitly include the DMF molecules into the calculations, modelling the liquid phase. Since these calculations require a significant amount of CPU time a compromise must be made as to the computational model chosen: the PBE+D3 functional has been selected with the DZVP basis set. As the energetics are particularly affected by the functional, additional static, gas-phase calculations with the M06 functional have been done to separate the effect of the functional from the effect of solvation. To this end five calculation models are compared: the M06 and PBE+D3 functionals are employed in pure gas-phase calculations (M06/6-31G* and PBE+D3/6-31G*, respectively) and also in calculations with the implicit solvent model (M06/6-31G*/i-solv. and PBE+D3/6-31G*/i-solv.). Note that the basis set is not varied. The results with larger basis sets can be found in ref. 10 and further discussion is given in the ESI.†

Table 1 Activation free barriers (kcal mol⁻¹) calculated with gas-phase calculations, implicit (i-solv., gas-phase) and explicit (e-solv., QMMM) solvent models

Reaction	(1)	(2)	(3)
M06/6-31G*	27.6	7.5	19.3
PBE+D3/6-31G*	23.2	4.3	11.8
M06/6-31G*/i-solv.	25.0	11.2	19.1
PBE+D3/6-31G*/i-solv.	20.4	5.5	11.2
PBE+D3/DZVP/e-solv.	26.3	6.4	11.5

---

Table 2 Reaction free energies (kcal mol⁻¹) calculated with gas-phase calculations, implicit (i-solv., gas-phase) and explicit (e-solv., QM/MM) solvent models

Reaction	(1)	(2)	(3)
M06/6-31G*	20.9	−22.6	15.6
PBE+D3/6-31G*	19.3	−26.6	8.2
M06/6-31G*/i-solv.	11.5	−15.2	15.6
PBE+D3/6-31G*/i-solv.	13.0	−21.7	6.9
PBE+D3/DZVP/e-solv.	24.8	−17.5	8.6

---

The first observation is that switching from the M06 to PBE functional provides extra stabilization to the activation energies both with and without solvent. This is a well-known effect of GGA functionals to stabilize delocalized situations (such as elongated bonds in transition states) due to the improper treatment of the electron self-interaction energy. The 4–8 kcal mol⁻¹ decrease in the barriers is therefore attributed to this effect. On the other hand, the functional switch does not have a uniform effect on the reaction free energies: without solvent the PBE functional yields extra stabilization, whereas in the presence of solvent the effect is not systematic. However, we can see that the trends among the models are not affected by changing the M06 functional to PBE: reaction 2 is the most favorable whereas reaction 1 is the least favorable both kinetically and thermodynamically.

Comparison of rows 1 and 3, as well as rows 2 and 4, in Tables 1 and 2 shows that applying the SMD model does not yield systematic changes. Indeed, the different initial and final charge distributions imply that the solvent effects are also different along the paths leading to the observed variations. Note however that for the SMD model the average deviation between calculated and measured solvation energies is around 4 kcal mol⁻¹; hence some discrepancies between the results obtained from different solvation models are expected.³⁰

Next we turn to the results obtained from QM/MM calculations. The free energy curves obtained from the simulations are shown in Fig. 1. We can see that for reactions 2 and 3 the results agree nicely with those obtained by the same functional (PBE+D3) and SMD model. On the other hand, the calculations yield a more destabilized TS and product state for reaction 1. The values obtained by the explicit solvent model are closer to those obtained without solvent. This variation partly can be attributed to the lack of diffuse functions which could otherwise stabilize the pure negative states. It is clear that the effect depends on how the environment is treated (no solvent, implicit or explicit). Again, we stress that the trends characterizing the paths are not affected by this difference.

Fig. 1 Free energy profiles obtained from the QM/MM thermodynamic integration calculations. The continuous lines are smoothing cubic splines for display purposes.

An important issue is how the structural parameters are affected by the solvent models. As an example, Table 3 compares the distance of the C–O bond closing the furan ring at the transition state. It is seen that the variations in the CO distance for a given reaction can be attributed to both the functional and the nature of solvation: switching from the M06 to PBE functional elongates the CO bond at the transition state in the presence of implicit solvent (rows 3 and 4 in Table 3) while changing from implicit to explicit solvation shortens the TS distances.

Table 3 Comparison of the forming C–O bond distances (in Å) at the TS obtained with gas-phase calculations, implicit (i-solv., gas-phase) and explicit (e-solv., QM/MM) solvent models

Reaction	(1)	(2)	(3)
a Two decimal digits are given because the TS position is determined from the maximum of the free energy curve obtained from the TI sampling which yields less accurate values.
M06/6-31G*	1.839	2.330	1.850
PBE+D3/6-31G*	1.831	2.288	1.903
M06/6-31G*/i-solv.	1.919	2.203	1.867
PBE+D3/6-31G*/i-solv.	1.961	2.392	1.958
PBE+D3/DZVP/e-solv.^a	1.79	2.31	1.90

---

Further insight into the nature of solvation with explicit DMF molecules can be obtained by analyzing the trajectories from the MD simulations. First, we note that the selected solvent (DMF) is a strongly polar molecule; therefore a considerable solute–solvent interaction is expected. In fact gas-phase pairwise interaction energies (Kohn–Sham electronic energies) are quite large: they span an interval from −16.4 kcal mol⁻¹ to −36.9 kcal mol⁻¹, with the actual values depending on the charge state of the solute (weaker interaction with the negatively charged solute, strong interaction with the positively charged solute and intermediate values for the neutral solute; for actual values see ESI†). To find specific interaction patterns between solute and solvent, the 2D radial and 3D spatial distribution functions (SDF) of the solvent DMF molecules around the solute reacting molecule have been calculated at each value of the reaction coordinate. Fig. 2 shows a representative radial distribution function of the solvent oxygen around the carbon atom bonded to silver (see Scheme 1) for reaction 2 at the TS ensemble. It is seen that the density fluctuations are small, and they essentially indicate the very limited structuring effect of the solute–solvent interactions. Similar observations can be made for the other stages of the reactions.

Fig. 2 Radial distribution function of the DMF molecules (represented by the O atom) around the solute reactant (represented by the carbon atom which remains bonded to the Ag ion along the reaction path) for reaction 2, calculated for the TS ensemble.

Fig. 3 displays the evolution of the minimal Ag–solvent distances at different stages of reaction 2. The three panels show large solvent mobility for the initial, TS and product stages: the fluctuation of the minimal Ag–O distance can be on the order of 0.5 Å. This range is due to the combined translational–rotational movements of the solvent molecules as inspection of the trajectories indicates. The panels indicate that the solvent molecules are not immobilized around the solute molecule but instead are collectively fluctuating around it. This is an important observation as it does not follow directly from the strong pairwise interactions occurring in gas phase (see earlier). We have obtained similar patterns for reaction 3 (see ESI†). Note that as reaction 1 does not feature Ag this analysis was not done in that case.

Fig. 3 Minimal Ag–solvent distances represented by the instantaneous d(Ag–O) values as a function of time during a 3 ps interval at different stages of TI simulation of reaction 2. Top panel: initial stage; middle panel: TS; bottom panel: product stage.

The 2D representation of the solvent distribution averages out a lot of information. In contrast, a SDF shows also the angular distribution of the time-averaged density of the solvent and can more effectively reveal regions around the solute molecule where the interaction between solute and solvent may induce structural changes, such as formation of stable solvation shells. Fig. 4 displays the calculated SDF for the TS ensemble of reaction 2 in two subsequent periods of 2.5 ps from two different viewing orientations, where the colours indicate the different time intervals. It is seen that no steady regions around the reacting molecule can be identified, because the green and violet regions (the averaged higher density regions in the subsequent periods) are hardly overlapping. Similar plots have been obtained for the ensembles of the other values of the reaction coordinate. This indicates that the observed higher density regions are the results of temporal solvent fluctuations which are weakly affected by the solute. In particular, the lack of stable regions around the different functional groups of the solute molecule (carbonyl, ester, aromatic rings) indicates that all the interactions are of electrostatic nature and no directional interactions such as H-bond can form between solvent and solute. This is not surprising as neither the solvent nor the solute molecules contain suitable H-bond donor groups. Note however that the large dielectric constant of the bulk DMF (36.7 at 298 K) and the large dipole moment (3.9 D) of the individual solvent molecules hint at considerable electrostatic interactions between solvent and the polarized solute molecules such as in the present case. Clearly however, while a quite compact and close solvation shell can be recognized around the solute, it is still quite mobile.

Fig. 4 Spatial distribution function of the central N atom of the solvent DMF molecules around the silver cation as the origin of the solute molecule calculated at the TS ensemble of the cyclization reaction 2 in two consecutive periods of 2.5 ps. The green and magenta isosurfaces represent ten times the uniform densities in the two time intervals. The two plots represent the same distribution seen from different angles.

Essentially the same conclusions can be drawn by inspecting the pseudo spatial distribution function (PSDF) of the solvent around the solvent molecule. In Fig. 5 the PSDF is plotted for the TS ensemble in two successive periods of 2.5 ps. To obtain the PSDF plots the solvent distribution has been cylindrically averaged out along the selected trajectory segment for the solute orientation shown in Fig. 2. A similar behaviour can be obtained for the other stages of the reactions. The plots show that no specific, stable solvation region is formed, as the excess densities are typical for the whole solvent region and these denser regions do not survive the selected short, 2.5 ps periods.

Fig. 5 Pseudo spatial distribution function of the central N atom of the solvent DMF molecules around the axis defined by the vertical C–C bond indicated on the plots calculated at the TS ensemble of the cyclization reaction 2 in two subsequent periods of 2.5 ps.

The analysis of the distributions of the solvent molecules indicates that the solvent is only a spectator and does not form any particular solvation pattern, i.e. the movements of the solvent molecules are not confined to specific regions around the solute. In other words, the movements of the solvent molecules can be averaged out and their individual, time-dependent effects can be replaced by ensemble-averaged effects. We thus see that the explicit inclusion of the solvent molecules does not bring any benefit for the calculations: the energetics and structural properties are very similar for both the explicit and implicit models and the solvent behaviour does not feature any direct mechanistic participation. This implies that the implicit solvent model is perfectly adequate in this case but for a fraction of the computational costs.

There are many examples where explicitly including the solvent into the calculation model has been essential to obtain appropriate chemical insight into reactivity issues.³¹ In these cases the solvent participates directly in the reaction by opening a new route for the reaction steps. Such situations are when the solvent has long-lived interactions with the reactants or it is itself a reactant,^7,31c,e or when it stabilizes/destabilizes selected states of the reaction (reactant state destabilization or TS stabilization). In the present case no direct solvent participation in the reaction occurs. The lack of a steady patterned solvation shell around the reacting solute implies that the effect and complexity of the solvent can be described by averaging out the distribution of the solvent molecules, which in turn is proved to be successfully approximated by implicit solvent models.⁵ This is the reason why the implicit solvation model works efficiently for this reaction step.

The above discussion provokes the question whether it can be predicted in advance if implicit solvation is sufficient in a given simulation. Although the present study discusses only three reactions with different charge states, a generalization is worth noting because this can save an enormous CPU time and cost: one can anticipate that implicit solvation works sufficiently in similar processes when no active solvent participation can be assumed in the reaction and no specific solute–solvent interaction is expected. This does not exclude considerable interaction between individual solvent–solute (as in the present case) but without inducing any kind of ordering in the solvation environment. It is important to see that the present TI simulation with explicit solvation for obtaining free energy barriers and reaction free energies tacitly assumes that the solvent included explicitly in the calculations equilibrates fully around the intramolecularly reacting solute molecule during the reaction step; this assumption may not always be true, e.g. for reactive radical intermediates.³² However this assumption is also an important ingredient of the implicit solvation models, i.e. it is a necessary condition both for implicit solvent models and for explicit models operating with equilibrium simulations.

Conclusions

Comparison of calculations with explicit and implicit solvent models on non-catalysed and Ag-catalysed ring-closing reaction steps featuring three different charge states leads to the following conclusions: (i) the effect of the solvent DMF (non H-bond forming solvent) can be adequately described by implicit solvent models; (ii) the simulations have demonstrated that the solvent does not form a stable solvation shell; (iii) site-specific interactions between the solvent molecules and the reactive solute do not play a role in the reaction paths investigated; and (iv) for the present reaction routes inclusion of the solvent effects does not change the reactivity trend obtained even without considering any solvent effect. The selected set of reactions represents several analogous situations in which metal cations catalyse organic coupling reaction steps. An important message of the calculations is that although explicit consideration of the solvent offers improvement in the molecular description, the simplified chemical model of the continuum solvation may be qualitatively correct even in the case of more polar solvents. Whenever there is an indication that the solvent does not participate explicitly in the mechanism and specific ordering of the solvent by site-specific interactions is not expected the implicit solvation model is a sensible choice. We note, however, that the above conditions are only necessary but not sufficient conditions for the proper choice of the solvation model because actual situations do not necessarily comply with other assumptions built in the implicit models. Still, the QM/MM approach is unavoidable if these two factors are expected or when comparison with experiments disagrees with the calculated mechanism.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The work was supported by NKFIH grant K-116034. Computational resources provided by the NIIF Supercomputer Center are acknowledged. We are grateful to the anonymous reviewers to point out important issues in the original version of the manuscript.

Notes and references

1. G.-J. Cheng, X. Zhang, L. W. Chung, L. Xu and Y.-D. Wu, J. Am. Chem. Soc., 2015, 137, 1706–1725.
2. T. Sperger, I. A. Sanhueza and F. Schoenebeck, Acc. Chem. Res., 2016, 49, 1311–1319.
3. M. Besora, P. Vidossich, A. Lledós, G. Ujaque and F. Maseras, J. Phys. Chem. A, 2018, 122, 1392–1399.
4. B. Mennucci, J. Phys. Chem. Lett., 2010, 1, 1666–1674.
5. J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3094.
6. J. Gao, Acc. Chem. Res., 1996, 29, 298–305.
7. A. Stirling, N. N. Nair, A. Lledós and G. Ujaque, Chem. Soc. Rev., 2014, 43, 4940–4952.
8. M. M. Heravi, M. Ghavidel and L. Mohammadkhani, RSC Adv., 2018, 8, 27832–27862.
9. C. He, S. Guo, J. Ke, J. Hao, H. Xu, H. Chen and A. Lei, J. Am. Chem. Soc., 2012, 134, 5766–5769.
10. J. Daru, Z. Benda, Á. Póti, Z. Novák and A. Stirling, Chem. – Eur. J., 2014, 20, 15395–15400.
11. (a) G. Fang and X. Bi, Chem. Soc. Rev., 2015, 44, 8124–8173; (b) G. Fang, X. Cong, G. Zanoni, Q. Liu and X. Bi, Adv. Synth. Catal., 2017, 359, 1422–1502.
12. (a) M. Chioua, E. Soriano, L. Infantes, M. L. Jimeno, J. Marco-Contelles and A. Samadi, Eur. J. Org. Chem., 2012, 35–39; (b) C. Liu, Y. Luo, W. Zhang, J. Qu and X. Lu, Organometallics, 2014, 33, 2984–2989; (c) E. Boz and N. Ş. Tüzün, Dalton Trans., 2016, 45, 5752–5764; (d) J. Ueda, S. Harada, H. Nakayama and T. Nemoto, Org. Biomol. Chem., 2018, 16, 4675–4682; (e) R. Roohzadeh, B. Nasiri, A. Chipman, B. F. Yates and A. Ariafard, Organometallics, 2018, 38, 256–267.
13. C. Reichardt and T. Welton, Solvents and Solvent Effects in Organic Chemistry, 4th edn, Wiley-VCH, 2011.
14. (a) The CP2K Developers Group, available at: http://www.cp2k.org; (b) J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128; (c) J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2013, 4, 15–25.
15. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
16. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
17. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
18. C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641–3662.
19. A. D. MacKerell, J. Wiorkiewicz-Kuczera and M. Karplus, J. Am. Chem. Soc., 1995, 117, 11946–11975.
20. T. Laino, F. Mohamed, A. Laio and M. Parrinello, J. Chem. Theory Comput., 2005, 1, 1176–1184.
21. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 14101.
22. (a) E. A. Carter, G. Ciccotti, J. T. Hynes and R. Kapral, Chem. Phys. Lett., 1989, 156, 472–477; (b) C. Chipot and A. Pohorille, Free Energy Calculations, Springer, Berlin, 2007.
23. For details see ESI†.
24. M. Brehm and B. Kirchner, J. Chem. Inf. Model., 2011, 51, 2007–2023.
25. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396.
26. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision E.01, Gaussian, Inc., Wallingford CT, 2013.
27. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2007, 120, 215–241.
28. Y.-F. Yang, G.-J. Cheng, P. Liu, D. Leow, T.-Y. Sun, P. Chen, X. Zhang, J.-Q. Yu, Y.-D. Wu and K. N. Houk, J. Am. Chem. Soc., 2013, 136, 344–355.
29. P. Xiao, H. Yuan, J. Liu, Y. Zheng, X. Bi and J. Zhang, ACS Catal., 2015, 5, 6177–6184.
30. (a) W. Wu and J. Kieffer, J. Chem. Theory Comput., 2018, 15, 371–381; (b) A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Chem. Theory Comput., 2009, 5, 2447–2464.
31. Selected examples: (a) J. Gao, Acc. Chem. Res., 1996, 29, 298–305; (b) B. K. Carpenter, J. N. Harvey and A. J. Orr-Ewing, J. Am. Chem. Soc., 2016, 138, 4695–4705; (c) A. Stirling, J. Phys. Chem. B, 2011, 115, 14683–14687; (d) S. Ahmadi, L. Barrios Herrera, M. Chehelamirani, J. Hostaš, S. Jalife and D. R. Salahub, Int. J. Quantum Chem., 2018, 118, e25558; (e) J. J. Varghese and S. H. Mushrif, React. Chem. Eng., 2019, 4, 165–206.
32. (a) B. K. Carpenter, J. N. Harvey and A. J. Orr-Ewing, J. Am. Chem. Soc., 2016, 138, 4695–4705; (b) A. J. Orr-Ewing, Chem. Soc. Rev., 2017, 46, 7597–7614.

---

Footnote

† Electronic supplementary information (ESI) available: Further comparisons of the methods, verification of the initial equilibria, details of the statistical error calculations, Cartesian coordinates of the stationary structures. See DOI: 10.1039/c9nj04003j

---