Role of Microsolvation and Quantum Effects in the Accurate Prediction of Kinetic Isotope Effects: The Case of Hydrogen Atom Abstraction in Ethanol by Atomic Hydrogen in Aqueous Solution

Role of Microsolvation and Quantum Effects in the Accurate Prediction of Kinetic Isotope Effects: The Case of Hydrogen Atom Abstraction in Ethanol by Atomic Hydrogen in Aqueous Solution Suraj Kannath, PawełAdamczyk, David Ferro-Costas, Antonio Fernandez-Ramos, Dan Thomas Major, and Agnieszka Dybala-Defratyka* Cite This: J. Chem. Theory Comput. 2020, 16, 847−859 Read Online ACCESS Metrics & More Article Recommendations * sıSupporting Information ABSTRACT: Hydrogen abstraction from ethanol by atomic hydrogen in aqueous solution is studied using two theoretical approaches: the multipath variational transition state theory (MPVTST) and a path-integral formalism in combination with free-energy perturbation and umbrella sampling (PI-FEP/UM). The performance of the models is compared to experimental values of H kinetic isotope effects (KIE). Solvation models used in this study ranged from purely implicit, via mixed−microsolvation treated quantum mechanically via the density functional theory (DFT) to fully explicit representation of the solvent, which was incorporated using a combined quantum mechanical-molecular mechanical (QM/MM) potential. The effects of the transition state conformation and the position of microsolvating water molecules interacting with the solute on the KIE are discussed. The KIEs are in good agreement with experiment when MP-VTST is used together with a model that includes microsolvation of the polar part of ethanol by five or six water molecules, emphasizing the importance of explicit solvation in KIE calculations. Both, MP-VTST and PI-FEP/UM enable detailed characterization of nuclear quantum effects accompanying the hydrogen atom transfer reaction in aqueous solution. 1. INTRODUCTION Hydrogen abstraction from ethanol by atomic hydrogen is a well-known reaction which is one of the most important steps in ethanol decomposition. 1−5 This type of hydrogen abstraction reaction is used as competition kinetic standards to determine relative reaction rate constants for different solutes where hydrogen is not released during the reactions. 2−5 Depending on temperature, this reaction can proceed via three different channels arising from the hydrogen atom abstracted from different positions within the ethanol molecule (methyl or methylene or the hydroxyl group). 6,7 However, at room temperature the reaction involving the abstraction from the closest carbon atom to the hydroxyl group (α-carbon) is much faster than the H abstraction from the hydroxyl or methyl groups. Therefore, only this channel needs to be considered for the reaction occurring in aqueous solution (reaction R1). CH CH OH H CH CHOH H 32 3 2 +→+(R1) Ethanol presents three distinguishable minima with the methyl group in various positions with respect to the hydrogen atom of the hydroxyl group. In two of the structures the methyl group is in the gauche configuration (g+ and g−), whereas in the other one it is in the alternate or trans (t) configuration. These three structures interconvert easily by internal rotations about the C−O bond. 8 One would expect also three transition state structures, but there are only two, with the methyl group in g (g+, or g−), depending on the hydrogen abstracted and in t configurations with respect to the OH group (Figure 1). Since the hydrogen abstraction generates an asymmetric carbon at the transition state in the αposition, there are two additional transition states when the hydrogen is abstracted from the other hydrogen atom of the α-carbon. The two transition states resulting from the H abstraction of one of the hydrogen atoms are configurational enantiomers of the transition states obtained from the abstraction of the other hydrogen. Therefore, the total rate constant for the process is twice the one obtained from the one that only considers the abstraction of one hydrogen atom (Figure 1). In chemical and biological reactions, kinetic isotope effects (KIEs) act as an important tool for explaining reaction mechanisms. 9−15 Comparison of KIEs originating from measurements with those determined computationally plays Received: August 2, 2019 Published: January 6, 2020 Articlepubs.acs.org/JCTC © 2020 American Chemical Society 847 https://dx.doi.org/10.1021/acs.jctc.9b00774 J. Chem. Theory Comput. 2020, 16, 847−859 This is an open access article published under a Creative Commons Attribution (CC-BY) License, which permits unrestricted use, distribution and reproduction in any medium, provided the author and source are cited. Downloaded via UNIV DE SANTIAGO DE COMPOSTELA on January 26, 2021 at 12:27:04 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. a crucial role in understanding, validating, or rejecting mechanisms of chemical or biological reactions. 16−19 KIEs are defined as the ratio of rate constants for two reactions which only differ in their isotopic composition. Three different isotopic substitutions along with the contribution of each conformer to the total rate constant and accompanying isotopic fractionation were studied (reactions R2,R3, and R4). CD CD OH H CD CDOH HD 32 3 +→+(R2) CH CH OD D CH CHOD DH 32 3 +→+(R3) CD CD OD D CD CDOD D 32 3 2 +→+(R4) Among the existing theories for calculating biomolecular rates, the transition state theory (TST) is one of the most widely used. 9,20,21 The main advantage of the TST over other theories is that it is easy to use and most suitable for analogizing the trends of a chemical reaction in terms of quantities which can be easily interpreted. 20,21 In the past few decades, there has been extensive research going on for exploring the parameters that control the chemical reaction rates. 1,22−25 Hydrogen abstraction reaction is a subject of intense theoretical research focused on seeking the most adequate treatment of quantum effects which constitute the origin of KIEs and are necessary to describe the behavior of light elements like hydrogen and its isotopes. One of the theoretical approaches is the variational transition state theory (VTST) 22,24,26,27 and its recent developments, such as, e.g., multipath VTST (MPVTST), 8,28−30 which allows for taking care of reactions with multiple conformers. This methodology also enables the incorporation of zero-point energy and multidimensional tunneling variational effects as well as multiple conformations in the determination of thermal rate constants. All these corrections can be used to make more accurate predictions of KIEs. 1,8 It is well-known that zero-point energy and nuclear quantum mechanical effects are essential when studying H transfer reactions. 32,33 When dealing with solvent effects it is more common to model them using a continuum solvent methodology, so explicit solvation is typically not investigated within VTST calculations 1,8,31 (however, the so-called ensemble averaged VTST with multidimensional tunneling, EA-VTST/MT, which allows for incorporation of explicit environment, was regularly used a few years ago to study numerous enzymatic reactions 32−43 ). This implicit solvent assumption may lead to severe simplification, in particular in cases where specific solute−solvent interactions play an important role. Free radical kinetics in solution can be a good example of such cases as even more sophisticated treatment of solvation (by including equilibrium and nonequilibrium effects) is not always sufficient for reproducing all kinetic isotope effects measured on a studied reaction. 44,45 Among available methods which conveniently allow for the presence of explicit solvent molecules in bulk we consider Feynman path integrals (PI). In particular, using a combined quantum mechanical and molecular mechanical (QM/MM) potential in conjunction with PI is a powerful tool for determining KIEs. 46,47 The integrated discretized PI method and free-energy perturbation/umbrella sampling (PI-FEP/ UM) approach in molecular dynamics simulations (MD) have been found useful for a variety of applications including H transfer reaction in solution and biological systems. 48−52 It is important to notice that both chosen theories, MP-VTST and PI-FEP/UM, have their limitations. For instance, PI-FEP/UM calculations provide total nuclear quantum effect (NQE) but do not easily allow dissection of tunneling contributions, while MP-VTST calculations naturally provide tunneling contributions. Within the PI-FEP calculations it is straightforward to take into account a large number of solvent molecules (treated classically), while in MP-VTST specific solvent molecules must be manually placed around the solute, possibly within a continuum solvent medium. This makes both theories complementary and allows this approach to provide a deeper description of the reaction under study. In this work, we used both VTST and PI-FEP to predict hydrogen KIEs on H atom abstraction from an ethanol molecule in aqueous solution. A very recent computational study on this reaction by some of us, using a continuum model of solvation, 8 provided results very close to the experimental values in the case of R2 substitution. However, in the case of R3, no isotope effect as opposed to an inverse kinetic isotope effect was predicted, and R4 substitution resulted in values slightly lower than experimental ones. Hence, these earlier results leave room for improvement, 7,8 for instance, by including explicit solvation effects. The goal of this work is to evaluate and understand the influence of explicit solvation in the calculation of KIEs in the H-abstraction in ethanol. Using MP-VTST we investigate the number of specific microsolvating water molecules required, within a continuum environment, to reach converged KIE values. These values are compared with those obtained from PI-FEP/UM simulations, which include fully explicit solvent treatment. Such a study using VTST and PI allows a careful comparison of two widely used complementary theories and the role of solvation. 2. THEORY Path-Integral Free Energy Perturbation/Umbrella Sampling (PI-FEP/UM). In a system where the entire environment is treated explicitly the rate constant of reaction can be determined using the path-integral quantum transition state theory (QTST) 53,54 Figure 1. Ethanol (top panel) and transition state (bottom panel) conformations. H[1] and H[2] indicate the hydrogen atom being abstracted and the abstracting hydrogen atom, respectively. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://dx.doi.org/10.1021/acs.jctc.9b00774 J. Chem. Theory Comput. 2020, 16, 847−859 848 khe 1F QTST qm β =β−Δ‡ (1) where his Planck’s constant, β=1/kBTwhere kBis Boltzmann’s constant, and Tis the temperature. ΔFqm ‡is the free energy of activation which can be obtained from the quantum mechanical potential of mean force (PMF, wqm(z)) Fwzwz() ( ) qm qm qm R Δ =− ‡‡ (2) where zis the centroid reaction coordinate, and z‡specifies its value at which wqm(z) has its maximum value, and zRdenotes the coordinate at the reactant state. Feynman path-integral simulations allow for determining wqm(z), on the one hand, and incorporating nuclear quantum effects, on the other hand. 55 Many practical methods have been developed over the years to determine the quantum mechanical reaction rate constant. Those procedures include centroid molecular dynamics 56 or ring-polymer molecular dynamics. 57 The methodology that is used in this work is an alternative approach within which the quantum mechanical PMF (wqm(z)) is determined through a double averaging procedure and is called quantized classical path (QCP). 58,59,72−75 Using this methodology KIE can be calculated as the ratio of two rate constants 58 k k e e K IE wz wz wz wz L H () ( ) () () qm LLqm LL R qm HHqm HH R == β β −{ −} −{ −} ‡ ‡(3) where L and H denote light and heavy isotopologues. According to this equation, KIE is predicted based on the PMFs computed separately for each isotope. However, this approach has been shown to provide large statistical fluctuations leading to inaccurate results. By expressing KIE in terms of the partial partition functions, the Qqm ratio, which can be determined directly through the free energy perturbation (FEP) theory by perturbing the mass from one isotope into the other, one can obtain KIE from a single pathintegral simulation: 58 k k Qz Qz Qz Qz K IE () () () () L H qm LL qm HH qm HH R qm LL R ==     ‡ ‡ Ä Ç Å Å Å Å Å Å Å Å Å Å Å Å É Ö Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ä Ç Å Å Å Å Å Å Å Å Å Å Å Å É Ö Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ(4) Since no modifications to the original methodology are provided in this work, the more interested reader is referred to the original papers as well as two recent articles in which the entire approach has been reviewed. 46,60 The FEP theory was also applied by others, 61 for instance in the SC-FEP scheme by Markland and Ceriotti. 62 Approaches based on FEP are successors of related alchemical transformation - thermodynamic integration (TI). 63,64 TI has problems of its own like integration error coming from the fact that the mass is discretized. One of the very recent approaches allowing to either reduce or remove that error was introduced by Karandashev and Vanicek 65 to accelerate the calculation of equilibrium isotope effects. Multi-Path Variational Transition State Theory (MPVTST). The canonical version (CVT) 66 of MP-VTST including small-curvature (SCT) multidimensional corrections for tunneling 67 can be written as kT TkT() () () MP CVT/SCT MP CVT/SCT MS TST γ=⟨ ⟩ ‐‐‐ (5) where kMS‑TST(T) is the multistructural transition state theory rate constant, which is given by kTBT QT QT () () () () MS TST MS HO, Et MS HO = ‐‐‡ ‐(6) where QMS‑HO,‡(T) and QEt MS‑HO(T) are the multistructural rovibrational rigid-rotor harmonic-oscillator partition functions that incorporate all transition state and ethanol minima (conformations), respectively; B(T) includes the electronic and translational partition functions for reactants and the transition state, respectively. Notice that for the cases in which ethanol is microsolvated the reactants and transition state partition functions include also the surrounding water molecules. In the case of ethanol, microsolvation affects the stability of the three minima of ethanol, because the distribution of the water molecules makes the g+ and g− configurations slightly unequal. The factor ⟨γMP‑CVT/SCT(T)⟩is an average of the canonical variational (Γ) and small-curvature tunneling (κ) corrections calculated with information obtained from each of the reaction paths that start at each of the transition state conformers. For a reaction with n‡transition states, which are conformational isomers TTQ T QT () () () () i n ii MP CVT/SCT 1 CVT/SCT RR HO, MS HO, γγ ⟨ ⟩= ∑ −= ‐‡ ‐‡ ‡ (7) where Qi RR‑HO,‡is the rigid-rotor harmonic-oscillator partition function of the i-th transition state structure, Qrv MS‑HO,‡(T)is the multistructural rovibrational partition function, and TT T() () () iii CVT/SCT CVT CVT/SCT γ κ=Γ (8) The rate constant for the passage through a given transition state is kT TkT TkT TTkT () () () () () () () () iiii i iii CVT/SCT CVT/SCT TST CVT/SCT CVT CVT/SCT CVT TST γκ κ == =Γ (9) The variational Γi CVT(T) coefficient of eq 8 accounts for recrossing effects. It is given by the ratio between the canonical variational theory (CVT), 66 ki CVT(T), and the transition state theory (TST), ki TST(T), thermal rate constants TkT kT () () () i i i CVT CVT TST Γ= (10) The tunneling contribution κi CVT/SCT(T) was calculated using the small-curvature approximation. 67 Notice that ki CVT/SCT(T), ki CVT(T), and ki TST(T) include not only the partition function of the i-th transition state but also the MS-HO partition function of reactants. This is because the barrier heights for interconversion between conformers are much lower than the barrier height of reaction. Thus, the sum over ki CVT/SCT(T) is the MP-CVT/SCT rate constant, whereas the sum over ki TST(T) leads to the MS-TST rate constant. The total KIE, η, is given as the ratio between the MP-CVT/ SCT rate constants for the root (H) and the isotopically substituted species (D), i.e., by k k H MP CVT/SCT D MP CVT/SCT η = ‐ ‐(11) Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://dx.doi.org/10.1021/acs.jctc.9b00774 J. Chem. Theory Comput. 2020, 16, 847−859 849 where we have dropped the temperature dependence. The total KIE can be expressed as a sum of individual contributions of each transition state conformer: i n i 1 ∑ η η = = ‡ (12) where P iii ,D η η =(13) can be considered a weighted contribution of the i-th transition state to the KIE. The weighting factor is given by P Q Q i ii i n ii ,D ,D CVT/SCT ,D RR HO, 1,D CVT/SCT ,D RR HO, γ γ = ∑ ‐‡ = ‐‡ ‡ (14) and ηiis the individual transition state contribution to the KIE, which is given by the ratio of the individual CVT/SCT thermal rate constants k k i i i ,H CVT/SCT ,D CVT/SCT η = (15) The individual transition state KIEs can be decomposed into the translational, ηtrans, rovibrational, ηrv,i ‡, and variational and tunneling contributions, ηvtun,i: iiitrans rv, vtun, η ηηη=‡(16) The individual rovibrational contribution includes the MS-HO partition function of reactants Q Q Q Q i i i rv, ,H RR HO, ,D RR HO, Et,D MS HO Et,H MS HO η = ‡ ‐‡ ‐‡ ‐ ‐(17) whereas the variational and tunneling contributions are simply given by i i i vtun, ,H CVT/SCT ,D CVT/SCT η γ γ = (18) Anharmonic effects due to the hindered rotors can be also included, but we have encountered from a previous work 8 that in this case its contribution to the KIE is negligible. 3. COMPUTATIONAL DETAILS PI-FEP/UM Simulations. In the case of PI-FEP/UM simulations, the QM atoms (ethanol molecule and hydrogen atom) were modeled using the AM1 Hamiltonian 68 and the Minnesota DFT MPWB1K functional 69 along with the 6- 31+G(d,p) 70−74 split valence-ζbasis set. This functional has been recommended for thermochemistry and thermochemical kinetics 75,76 and has been shown to perform well for hydrogen transfer reactions. 28,77,78 However, in order to justify its use for the current system, we have performed some tests with other methods, and the results are presented in the SI. Among tested electronic structure methods two functionals MPWB1K and M05-2X 76 resulted in the closest agreement with the electronic barrier of 7.6 kcal mol−1derived from the experimental activation energies and the ab initio calculations reported by Lossack et al. 7 (Table S1). Furthermore, the MPWB1K functional provided KIEs for the bare model, using the MPCVT/SCT calculation described in detail below, that better reproduced the experimental magnitudes (Table S2). The solvent region was represented explicitly by 1080 TIP3P 79 water molecules. Stochastic boundary conditions (SBC) 80 were employed (Figure 2). The solvated system was heated up to 298 K for 20 ps and then equilibrated for 30 ps. The reaction coordinate (z) was defined as the difference between the C−H1bond cleavage and the H1−H2bond formation. Such a simple choice of the reaction coordinate has been shown to be equally effective in describing the reaction where hydrogen species is transferred between two heavy atom centers as other more sophisticated reaction coordinate definitions. 81,82 In order to obtain the Potential of Mean Force (PMF), umbrella sampling 83,84 simulations were used. A total of 17 regions (windows) were needed to cover the entire reaction range of interest. Each individual window was equilibrated for 100 ps, and data collection was performed for subsequent 200 ps. The trajectories resulting from umbrella sampling simulations were later introduced to path-integral calculations leading to direct evaluation of the path-integral quantities using the QCP approach. Isotopic substitution was introduced by a free energy perturbation method, 85 and KIEs were calculated using eq 4. In order to perform PI calculations, the three key atoms (namely the C, H1, and H2atoms) were quantized with various numbers of beads (P= 8, 16, 32, 64, and 128). Around 10 000 classical configurations were used to perform path-integral simulations for which each isotope was combined with 10 path-integral steps per classical step. To estimate the quantum contributions, an enhanced form of QCP 86−88 known as the bisection sampling (BQCP) 89,90 and staging sampling methods (SQCP) 91−93 were used. All PI simulations were carried out using the CHARMM software. 94 In the case of QM(DFT)/ MM simulations, CHARMM interfaced with Q-Chem 95 was used. MP-VTST Calculations. VTST calculations were performed using the same combination of the density functional and basis set as in the case of PI-FEP/UM QM(DFT)/MM calculations (i.e., MPWB1K/6-31+G(d,p)). In this work we enlarged the previously studied bare model by explicitly adding 1, 2, 5, and 6 water molecules to study the effect of solvation. To this end a potential energy surface (PES) scan was performed using the Gaussian G09 software 96 to locate the feasible positions of water molecules near the ethanol molecule. Solvation models containing 5 and 6 water Figure 2. QM/MM model of H abstraction from ethanol. Water molecules are shown in wire representation; the QM region is shown as a ball-and-stick model. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://dx.doi.org/10.1021/acs.jctc.9b00774 J. Chem. Theory Comput. 2020, 16, 847−859 850 molecules were constructed from clusters of 6 and 7 water molecules, respectively, taken from the Cambridge Cluster Database. 97 The IEFPCM continuum solvent model 98 along with its default options was also used to include additional solvent effects. All geometry optimizations were carried out using very tight convergence criteria and ultrafine grid integration. The Gaussian 09 software was used to obtain the optimized geometries, energies, and Hessians. Harmonic normal-mode frequencies were scaled by a factor of 0.964. 99 The minimum energy path (MEP) was followed starting at each of the n‡transition states with the objective of calculating the variational and tunneling effect contributions to the MPCVT/SCT thermal rate constant. The MEP was obtained using the Page-McIver algorithm 100 in mass-scaled coordinates (scaling mass of 1 amu), with a step size of 0.005 au, and Hessians were calculated every 10 steps. The frequencies along the reaction path were projected using redundant internal coordinates. 101 The Pilgrim 102 software package interfaced with Gaussian 09 was used for calculation of the multipath canonical variational transition state rate constants with small curvature multidimensional tunneling corrections (MP-CVT/ SCT). 103 4. RESULTS AND DISCUSSION Full Explicit Solvation Model. The classical potential of mean force obtained for the reaction between an ethanol molecule and a hydrogen atom using the AM1 Hamiltonian for the QM region along the defined reaction coordinate, z,is presented in Figure 3. The free energy barrier was around 3 kcal mol−1, which is much lower than the electronic barrier height of 7.6 kcal mol−1estimated based on the measurement of activation energy and theoretical prediction reported by Lossack et al. 7 TheresultedKIEsobtainedfromthe subsequent PI-FEP/UM calculations for a different number of beads and two different sampling schemes are shown in Table 1. Both sampling schemes give similar results. The trends of normal KIEs for reactions R2 and R4 were preserved in all the calculations and all methods that were used to calculate the KIEs. In the case of R3 substitution, almost all calculations provided normal KIEs which do not agree with the experimental data. However, the magnitudes for the R4 substitution were significantly smaller for AM1 than the experimental ones. Since the results have converged at 32 beads, using a higher number of beads did not seem necessary. Since the barrier height from the QM(AM1)/MM simulations, as well as the predicted KIEs, was not satisfactory, we also used a DFT level for the QM part. Hence, we performed again umbrella sampling simulations using QM(DFT)/MM. Each individual window was equilibrated for 5 ps, and data collection was carried out for the subsequent 27 ps. The free energy barrier height increased to 10 kcal mol−1as shown in Figure 3. The trajectories from each window were subsequently used in the PI quantities evaluation. To obtain quantum PMF in this case, only bisection sampling was used. However, the number of classical configurations used was limited to 2,500 due to the huge computational expense of the QM(DFT)/MM simulations and the additional cost of PI sampling (i.e., number of beads and MC steps). Based on the convergence behavior using AM1, we concludedthatitissufficienttouse32beadswith MPWB1K. The calculated KIEs using MPWB1K reproduced the experimental trends, although the deviation from experimental values is uneven. The best agreement was obtained for R2 (only deviates by 5%), a bit worse for R3 (22%), and the poorest for R4 (26%) (Table 2). Another important observation is that with the increasing number of beads, the predicted KIEs for R2 and R4 tend to increase and for R3 tend to decrease in most of the cases. Unfortunately, full convergence for this system is not attainable due to the high Figure 3. Classical potential of mean force obtained using the AM1 and MPWB1K methods to treat the QM region within the QM/MM simulations as a function of the reaction coordinate (z). Table 1. Hydrogen Kinetic Isotope Effects Obtained Using the PI-FEP/UM Method, Different Number of Beads, Different Sampling Algorithms, and the AM1 Hamiltonian for the QM Region Treatment PI-FEP/UM no. of beads R1/R2 R1/R3 R1/R4 bisection sampling 8 4.56 ±0.42 1.27 ±0.05 2.60 ±0.15 16 5.47 ±0.45 1.23 ±0.04 2.88 ±0.13 32 6.06 ±0.32 1.19 ±0.06 3.25 ±0.18 64 6.00 ±0.46 1.19 ±0.05 3.17 ±0.21 128 6.01 ±0.34 1.11 ±0.06 3.20 ±0.23 staging sampling 8 4.96 ±0.41 1.26 ±0.03 2.69 ±0.50 16 5.33 ±0.20 1.26 ±0.03 2.94 ±0.24 32 5.52 ±0.70 1.21 ±0.04 2.95 ±0.20 64 5.87 ±0.53 1.25 ±0.03 2.96 ±0.21 experimental a 7.38−0.85 +0.94 0.73−0.05 +0.06 6.80−0.98 +1.28 a Reference 7. Table 2. Hydrogen Kinetic Isotope Effects Obtained Using the PI-FEP/UM Method, a Different Number of Beads, and the MPWB1K Functional for the QM Region PI-FEP/UM no. of beads R1/R2 R1/R3 R1/R4 86.79 ±0.71 1.00 ±0.05 3.62 ±0.22 16 5.86 ±0.47 0.99 ±0.07 3.51 ±0.37 32 6.97 ±0.64 0.89 ±0.10 5.01 ±0.29 experimental a 7.38−0.85 +0.94 0.73−0.05 +0.06 6.80−0.98 +1.28 a Reference 7. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://dx.doi.org/10.1021/acs.jctc.9b00774 J. Chem. Theory Comput. 2020, 16, 847−859 851 cost of these simulations, although the variations in the results are within the statistical errors. The nuclear quantum mechanical free energy corrections such as the free energy difference and individual free energy contribution for different isotopic substitutions are provided in Figure 4. Classical and PI-FEP quantum free energy of activation for H abstraction from ethanol obtained using QM(MPWB1K)/MM simulations as a function of the reaction coordinate (z). Numbers indicate the values of free energies of activation. Figure 5. Snapshots of reactants (R), transition state (TS), and products (P) geometry resulting from the QM(MPWB1K)/MM model (only water molecules within a 4 Å distance from the ethanol molecule are shown for clarity). Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://dx.doi.org/10.1021/acs.jctc.9b00774 J. Chem. Theory Comput. 2020, 16, 847−859 852 Figures S1 and S2. Estimated free energy profiles for all studied reactions are shown in Figure 4. The decrease of free energy of activation (total quantum effects, ΔΔGqm ‡)is−2.5, −1.6, −2.2, and −1.8 kcal/mol in the case of reactions R1,R2,R3, and R4, respectively. Detailed structural analysis of specific interactions between the solute species and the surrounding water molecules at the transition state structures resulting from the simulations revealed that in the nearest neighborhood of the ethanol molecule, in particular its hydroxyl group, there are three water molecules out of which one is the hydroxyl proton acceptor and two (sometimes more) other water molecules donate their protons to the hydroxyl oxygen. Examples of such interaction patterns are shown in Figure 5. QM Microsolvation Models. Ethanol is an amphiphilic molecule with a polar or hydrophilic part corresponding to the OH group and a hydrophobic part involving the hydrocarbon chain. In principle, one would expect that if the reaction occurs in the hydrophobic part, the continuum model would account for most of the effects due to the solvent. However, our previous results that used this representation of the reaction were not completely satisfactory. In light of the PI calculations, which show a solvation shell around the hydroxyl group, some models were built to mimic the solvation of this polar group. Figure 5 shows that at least three water molecules are needed to microsolvate the hydroxyl group: two of them acting as hydrogen bond donors and one of them acting as a hydrogen bond acceptor. The smallest clusters with these characteristics, and at the same time including stabilization of the water cluster by additional hydrogen bonds, involve five (5W model) or six (6W model) water molecules. To analyze the effect of a single water acting as donor or acceptor, we also included models with one water, named 1W_d and 1W_a, in which the water molecule acts as a hydrogen bond donor and acceptor, respectively. Finally, we also added a model that was simply named as 2W, which is a combination of the 1W_d and 1W_a models. The transition state structures in the alternate configuration of the models and some of the reaction and intermolecular distances obtained at the MPWB1K/6-31+G(d,p) level are shown in Figure 6. Notice that the intermolecular distances between the water molecules and the OH group in the largest clusters are in good agreement with the PI-FEP calculations. The agreement is less satisfactory for the three atoms involved in the reaction. The QM(DFT)/MM model resulted in transition state geometries with C−H1and H1−H2distances which are too short and too large, respectively, when compared with the QM cluster calculations. The C−H1−H2bond angle is almost linear in the QM cluster (about 177 degrees) in all cases, whereas in the QM(DFT)/MM simulations this bond angle is reduced by almost 20 degrees. This may be the consequence of having all solvent molecules treated only classically in the QM(DFT)/MM model whereas in the microsolvated QM model they were computed at the same level of theory as the solute species or the result of thermal effects in the QM(DFT)/MM molecular dynamics simulations. Hereafter, the discussion is centered on the alternate configuration of the transition states (the one with the lowest energy), although the same arguments are also valid for the gauche configuration. The bond distances displayed in Figure 6 differ substantially in the 1W_a and 1W_d models. Both models show that the type of hydrogen bonding between the water molecule and the hydroxyl group has an important effect on the reaction site. The progress in the hydrogen abstraction coordinate at the transition state is greater in the 1W_d model. In fact, the reaction proceeds with lower barrier if the water molecule acts Figure 6. Transition state structures of the alternate configuration resulting from all microsolvation models and key interatomic distances (in Å). Values for the gauche configuration are shown in brackets. QM/MM distances are provided in the gray rectangular. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://dx.doi.org/10.1021/acs.jctc.9b00774 J. Chem. Theory Comput. 2020, 16, 847−859 853 as a hydrogen bond acceptor (6.37 vs 5.35 kcal·mol−1for the bare and 1W_a models, respectively). The opposite effect is observed when the water molecule acts as a hydrogen bond donor (6.37 vs 7.24 kcal·mol−1for the bare and 1W_d models, respectively, Table S4). The hydrogen bond distances in the 1W_a and 1W_d models for the trans reactants are 1.865 and 1.848 Å, respectively, whereas they are 1.800 and 1.900 Å at the corresponding transition states. Consequently, for the 1W_a model, the hydrogen bond is strengthened, which is in accord with the reduction of the reaction barrier. These two models and the 2W models are useful to look for trends but are inadequate to mimic the network of hydrogen bonds between the water molecules and the OH group. A good representation of this hydrogen bonding network can be achieved either with the 5W model or with the 6W model. With respect to the progress of the reaction coordinate at the transition state, these two models resemble the 1W_d model. In general, they present very similar geometric features, and the barrier heights for hydrogen transfer are almost the same (6.92 and 6.98 kcal·mol−1for the 5W and 6W models, respectively, Table S4). Therefore, it is expected to obtain similar thermal rate constants and KIEs for both clusters. Thermal rate constants and KIEs were calculated for all models using MP-CVT/SCT. Notice that the calculated MPVTST rate constants make use of the equilibrium solvation path (ESP) 72 approximation, as in ref 8. The variational and quantum effects due to the hydrogen abstraction reaction are plotted in Figure 7. The KIEs of Table 3 were factorized as indicated in eq 5. A full account of the results is given in the Supporting Information (Table S5). The calculated MP-CVT/SCT thermal rate constants (Table S3) are somewhat lower than the experimental values, being that the largest discrepancy between both sets is about a factor of 5. The agreement is not perfect, but in this work we concentrate on KIEs, which to a certain extent are independent of the barrier height (usually the main reason for disagreement between theory and experiment), rather than pursuing a perfect match between rate constants. An important issue that arises when calculating thermal rate constants using VTST is the choice of coordinates when evaluating certain quantities along the MEP. It has been pointed out that the frequencies calculated along the MEP in curvilinear internal coordinates provide more physical results than Cartesian coordinates. 69 This is because the projected frequencies influence the free energy profile, which is used to calculate the variational effects. Additionally, they are also used to evaluate the ground-state vibrationally adiabatic potential, which in turn is used to evaluate the SCT tunneling transmission coefficient. Therefore, in this work we used redundant internal coordinates. The effect of the choice of coordinates over the variational and tunneling effects can be seen in Tables S9 and S10 for the 5W model. Cartesian coordinates led to a huge underestimation of KIEs on the R2 and R4 substitutions, whereas they produced a more pronounced (more inverse) KIE for R3. The contribution of each of the transition state configurations to the total KIE ηiis given by eq 13, where the sum extends over the alternate and gauche configurations; Pi,D (eq 14) is the weight of each transition state, with the condition ∑iPi,D = 1. This weight is almost preserved for the bare model and the 5W and 6W models (especially between the bare and the 6W models). Therefore, discrepancies in the KIEs should be attributed to one or several terms in which the total KIE was factorized rather than to the contribution of each configuration. In this case, the slight increase in the R1/R2 and R1/R4 ratios in the clusters is mainly due to the increment of the rovibrational contribution to the KIE. The rovibrational contribution to the KIEs largely depends on the frequency differences between the C−H and C−D breaking bonds. At reactants, this difference translates into a large value for the ratio between partition functions (being that the one for the deuterated bond is the largest). This mass effect also influences the transition state partition functions in the opposite direction but to a lesser extent, so the rovibrational KIE is normal. In the case of R1/R3, the transferred atom is the same, and the aforementioned effect does not influence the partition function of reactants. However, the formation of a H−H bond vs the formation of a H−D bond at the transition states leads to an inverse rovibrational KIE. One of the main failures of the bare model was its inability to reproduce the observed inverse KIE for reaction R3, despite the strong inverse rovibrational contribution. However, both 5W and 6W models overcome this limitation. In this case, the discrepancy can be traced to the term that accounts for deviations from the transition state theory (variational effects Figure 7. Plots of the κCVT/SCT and ΓCVT coefficients calculated for the alternate configuration of each microsolvation model of each reaction. Table 3. Total MP-CVT/SCT KIEs Obtained at 298 K reaction model R1/R2 R1/R3 R1/R4 bare model 7.00 1.02 4.85 1W_d model 8.48 0.93 5.24 1W_a model 6.31 1.48 5.18 2W model 6.61 1.01 4.85 5W model 7.78 0.86 5.29 6W model 7.69 0.85 5.22 experimental a 7.38−0.85 +0.94 0.73−0.05 +0.06 6.80−0.98 +1.28 a Reference 7. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://dx.doi.org/10.1021/acs.jctc.9b00774 J. Chem. Theory Comput. 2020, 16, 847−859 854 and tunneling). An inspection of Figure 7 shows that the tunneling factor remains practically the same, but variational effects (recrossing) are less pronounced in the cluster models, leading to the inverse total KIE. These results indicate that to reproduce all the experimental KIEs it is necessary, on the one hand, to explicitly solvate the polar part of the molecule by forming a network of hydrogen bonds. An incomplete solvation (1W_a and 2W models) is not effective and leads to inferior results compared to the bare model. The 1W_d model, although still insufficient, provided better results. On the other hand, the methodology should incorporate recrossing and quantum effects, as is the case in MP-CVT/SCT. For instance, if the microsolvation model is incorporated within the conventional transition state theory framework, then the calculated KIEs on reactions R3 and R4 provide poor results, i.e., KIEs are overestimated in the former case and underestimated in the latter (Table S3). When applying MP-CVT/SCT, tunneling effects dominate over the recrossing, so reactions R1 and R3 are sped up with respect to R2 and R4, and the KIEs progress in the right direction (Tables S5, S6, and S7). Within the QM/MM simulations it is possible to estimate the overall energetic barrier reduction (so-called total quantum effect, ΔΔGqm ‡) caused by the quantum nature of either protium (H) or deuterium (D) transfer in all isotopic scenarios. The largest effect was observed for reaction R1 (2.5 kcal·mol−1), then for R3 (reduction by 2.2 kcal·mol−1), and less than 2 kcal·mol−1for R2 and R4 (1.6 and 1.8, respectively). As the extent of lowering the free energy of activation is not necessarily correlated with the magnitude of tunneling transmission coefficients, in VTST this decrease in value can be approximately estimated by the factor GG G ZP E qm MP CVT/SCT MS TST Δ Δ≈Δ −Δ +Δ ‡‐ ‐ (19) The first two terms of the rhs of eq 19 take into account the variation of free energy due to the variational and quantum effects (Table S11), and the last term includes the variation in the zero-point energy. The results are 2.35 and 2.78 kcal·mol−1 for R1 and R3, respectively, and 1.15 and 1.74 kcal·mol−1for R2 and R4, respectively. There is a qualitative agreement between the PI-FEP and VTST results, pointing toward the importance of quantum effects in this hydrogen abstraction reaction. This dual approach allowed for prediction of the KIEs for the reaction without losing much accuracy. By comparing H-KIEs obtained using 5W and 6W models and the MP-VTST approach with a full QM(DFT)/MM model and PI-FEP/UM approach, we conclude that both methods provide very good agreement with experiment, likely within the accuracy of the applied DFT method. Interesting alternatives to the MP-VTST and PI-FEP/UM calculations to account for microsolvation effects, usually performed at much lower temperatures than 298 K, are the ab initio molecular dynamics (AIMD) and the ab initio pathintegral (AIPI) simulations. 104 The latter treats the atomic nuclei as quantum degrees of freedom, and it has been, for instance, applied to the study of microsolvated HCl by water clusters 105 or to the microsolvation of protonated methane by molecular hydrogen. 106 In both research works, the authors stressed the importance of including quantum effects in the calculations, an aspect of the dynamics which is also of relevance to this study. However, dynamics quantum effects due to the solvent were not included in our study. This may be one of the reasons, together with the uncertainties in the electronic structure calculations, that prevents our models from providing a better match with experimental KIEs. MP-VTST is computationally cheaper than PI-FEP/UM for the current system, which has a small solute, which we surrounded by a few, selected water molecules, embedded in a continuum solvent. For reactions in which to locate all solvent configurations (as it may be the case in some enzyme or solution phase reactions) is computationally challenging, PI offers a good alternative to MP-VTST. However, in such cases semiempirical methods used to treat the QM region should rather be reparametrized 107−109 in order to describe the reaction under study at a higher accuracy. 5. CONCLUSIONS In the present study 2H KIEs were calculated for four different substitution scenarios using two theoretical approaches: a pathintegral formalism in the form of centroid path integral and free-energy perturbation-umbrella sampling (PI-FEP/UM) and the multipath variational transition state theory (MP-VTST). The former method enabled us to treat the solvent explicitly by using a QM/MM protocol and provided an educated guess for the construction of the water cluster models. It was observed that an incomplete solvation of the OH group did not improve the predictions compared to a continuum solvation approach. Rather, a complete solvation shell is required to reach converged KIE results. The latter method incorporated variational and quantum effects, such as tunneling, for multiple reaction paths in the evaluation of the KIEs. It was shown that both microsolvation and VTST corrections are needed in order to achieve a satisfactory comparison with experimental KIEs. Both PI-FEP/UM and microsolvation MP-VTST approaches provide results in good agreement with experiments. ■ASSOCIATED CONTENT * sıSupporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.9b00774. Comparison between several electronic structure methods for bare model; ratios of quantum mechanical partition functions; total free energy contributions for unsubstituted and substituted reactions; MS-TST, MPCVT, and MP-CVT/SCT thermal rate constants; classical barrier heights (V0 ‡) also with differences in zero-point energies included (ΔH0 ‡); variational and tunneling transmission coefficients; and key geometrical parameters for localized TSs (PDF) ■AUTHOR INFORMATION Corresponding Author Agnieszka Dybala-Defratyka −Institute of Applied Radiation Chemistry, Faculty of Chemistry, Lodz University of Technology, 90-924 Lodz, Poland; orcid.org/0000-0002-8939-2279; Email: [email protected] Authors Suraj Kannath −Institute of Applied Radiation Chemistry, Faculty of Chemistry, Lodz University of Technology, 90-924 Lodz, Poland Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://dx.doi.org/10.1021/acs.jctc.9b00774 J. Chem. Theory Comput. 2020, 16, 847−859 855