Reaction Dynamics of Cyanohydrins with ... - ACS Publications

Article Cite This: J. Phys. Chem. A XXXX, XXX, XXX−XXX

pubs.acs.org/JPCA

Reaction Dynamics of Cyanohydrins with Hydrosulfide in Water Stéphanie Valleau*,†,‡ and Todd J. Martínez*,†,‡ †

Department of Chemistry and The PULSE Institute, Stanford University, Stanford, California 94305, United States SLAC National Accelerator Laboratory, Menlo Park, California 94025, United States

‡

Downloaded via NOTTINGHAM TRENT UNIV on August 10, 2019 at 02:25:11 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: We studied the reaction dynamics of a proposed prebiotic reaction theoretically. The chemical process involves acetone cyanohydrin or formalcyanohydrin reacting with hydrosulfide in an aqueous environment. Rate constants and populations of reactant and product bimolecular geometric orientations for the reactions were obtained by using density functional theory for the energies, transitionstate theory for the rates, and matrix exponentiation as well as the hybrid tau-leaping algorithm for the population dynamics. The role of including the solvent explicitly versus implicitly was also investigated. We found that adding explicit water or hydrogen sulfide molecules lowers the activation energy barrier and leads to a more efficient reaction pathway. In particular, hydrogen sulfide was a better proton donor. Finally, we studied the role of including more than one reactant and product bimolecular orientation geometry in the dynamics. Including all initial pathways, reactant to reactant, product to product, reactant to product, and product to reactant led to a larger reaction rate constant compared to the single minimum energy pathway approach. Overall, we found that most reactions which involve formalcyanohydrin occur more rapidly or at the same speed as reactions which involve acetone cyanohydrin at room temperature.

■

INTRODUCTION Various reaction mechanisms for the prebiotic synthesis of nucleotides and amino acids have been studied and proposed by the scientific community.1−5 Recently, Sutherland et al.6 suggested a set of reactions which occur in the presence of hydrogen cyanide, HCN, and hydrogen sulfide, H2S, in a multiple pool geochemical scenario. The need for more than one pool comes from the different chemistry and reagents required in some steps of the process, for example the simultaneous presence of HCN and H2S impedes some reactions from occurring. The mechanism of many of the proposed reactions is unknown and it remains to be investigated whether they could have occurred subsequently and/or spontaneously with the periodic delivery of some reactants from flow chemistry. The optimal conditions, such as temperature, solvent, and external pressure, for each of the proposed reactions are also currently unknown. Nonetheless, these reactions provide a pathway which does not require the presence of any pre-existing external enzyme or catalystfor the generation of amino acids and nucleotides, stepping stones to the formation of RNA and DNA, the pillars of life in all organisms. In this context, a theoretical study of the kinetics can help provide information on the best conditions for the dynamics of each reaction. Here, we studied one specific reaction (Figure 1) from the network for the synthesis of amino acids,6 which occurs at room temperature in water. Our goal is that of understanding what level of theory is required for single reaction kinetics in order to eventually describe the kinetics of the overall network. We considered two options to include the solvent in the kinetics calculations, including hybrid implicit/explicit solvent molecules in the reactants with the implicit conductor-like screening model (COSMO) versus using only COSMO for the solvent.7 © XXXX American Chemical Society

Indeed, when protons are involved in the reaction (see e.g., refs8−10), the explicit inclusion of the solvent may lead to a more efficient pathway and a larger rate constant. We also investigated the influence of proton donors (H2S vs H2O) on rate constants. Much work on the theory and calculation of rate constants has been carried out in the last century.11,12 We recall the transition state theory (TST)13,14 and variational TST;15 both classical theories, mostly valid at high temperatures when tunneling is not important. At lower temperatures,16 approaches such as the quantum instanton17,18 or the calculation of the rate constant from the flux−flux or flux−side correlation functions19 take tunneling into account. Given that the reactions were carried out experimentally at 300 K and that the molecules are quite large, we opted to use TST for the rates with the Wigner correction to account in part for tunneling.20 TST is computationally less demanding given that the rate is estimated based on a transition state located on a single one-dimensional minimum energy path (MEP) connecting reactants to products. Therefore, TST avoids the cost of computing a large portion of the potential energy surface. Furthermore, one does not need to run any dynamics to obtain the rate with TST. By starting from the optimized geometry of reactants and products, we employed the simplified string approach21 as implemented in the TeraChem software package22−25 to find a MEP connecting reactants to products. Other approaches32 such as the nudged elastic band or the growing string method can also be used.26,27 With an optimized MEP and first order saddle point (transition state), it is then straightforward to compute the rate constant. Received: June 16, 2019 Revised: July 25, 2019 Published: July 26, 2019 A

DOI: 10.1021/acs.jpca.9b05735 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

mechanisms and rates in each case, and we also compared the roles of hydrogen sulfide and water as proton donors in the kinetics.

■

METHODS AND COMPUTATIONAL DETAILS With the goal of exploring the kinetics of the reactions in Figure 1 for multiple bimolecular reactants and product orientation geometries without having to compute a full potential energy surface, we chose to compute rates using TST with the Wigner correction for tunneling.13,14,20,31 The MEPs were obtained using the simplified string method21 as implemented in TeraChem.22−25 These rates were then combined to define a system of coupled ODEs for the overall reaction kinetics. The solution of this system of coupled ODEs was obtained using matrix exponentiation, VODE integration with backward differentiation formulas,28 and the hybrid SSA tau-leaping method.30,33 The overall procedure we followed is shown in the flow chart in Figure 2 and a detailed description of each part will be provided in the following subsections.

Figure 1. Reactions occurring in aqueous solvent studied in this work. Two reactants are considered, formalcyanohydrin (type I, R = H) and acetone cyanohydrin (type II, R = CH3). We consider the reaction of each of these two reactants with hydrosulfide (I and II) as well as their reaction with hydrosulfide and a water molecule (Iw and IIw), their reaction with hydrosulfide and a H2S molecule (Ihs and IIhs), and their reaction with hydrosulfide and both hydrogen sulfide and water (Iw+hs and IIw+hs). For reactions Iw, IIw, Ihs, IIhs, Iw+hs, and IIw+hs the product sulfur group is SH, whereas for the first two reactions, the S− anion remains unprotonated. The numbers of water or H2S molecules added in each reaction are indicated in the table above.

In general, when optimizing the geometry of multiple reactant species contemporarilyin our case, for example cyanohydrin and H2Svarious local minima can be identified. Therefore, using a single MEP for the kinetics requires very good a priori knowledge of the optimal geometry of reactants (respective orientation, distance etc.) and of the products, as well as knowing how representative that path is of the overall dynamics. This leads to the question: how significant is the role of reactant and product orientation geometries on the kinetics? To answer this question, for the first two reactions in Figure 1, we considered multiple paths originating from a set of bimolecular reactant orientation geometries and leading to a set of bimolecular product orientation geometries. For all other reactions (Iw, IIw, Ihs, IIhs, Iw+hs, IIw+hs), one or at most two MEPs were determined and the corresponding TST rate constants were calculated. The kinetics of the first two reactions was thus represented by a system of coupled ordinary differential equations (ODEs), which included information on all reactant to product, product to product, and reactant to reactant paths. Given that the system of ODEs was quite stiff, as the rate constant values ranged over 30 orders of magnitude, we considered three approaches to solve the coupled ODEs: exponentiation of the rate matrix, variable coefficient ODE (VODE) integration28 with backward differentiation formulas, and the stochastic hybrid SSA tau-leaping method.29,30Here SSA stands for stochastic simulation algorithm. We thus obtained populations for the bimolecular reactant and product geometries as a function of time for the reactions. Overall in this work, we computed and compared the rate constants for reactions with formalcyanohydrin or acetone cyanohydrin, we evaluated the role of including bimolecular reactant and product orientation geometries in the reaction kinetics, we studied the role of the solvent by computing rates with and without explicit solvent and determined the reaction

Figure 2. Flow chart which summarizes the procedure we followed to obtain the reaction rate constants and population dynamics. RMSD is the root mean square deviation, VODE stands for variable coefficient ODE solver, BDF stands for backward differentiation formulas, and SSA stands for stochastic simulation algorithm.

Initial Conditions. The geometries of formalcyanohydrin, acetone cyanohydrin, water, hydrogen sulfide, and hydrosulfide were optimized in TeraChem22−25 with density functional theory (DFT), the B3LYP functional, the 6-31G* basis set, the D3 dispersion correction,34 and the conductor like screening model COSMO (ε = 78.35). With these optimized geometries, the bimolecular reactant and bimolecular product orientation geometries for reactions I and II (see Figure 1) were obtained as following. We generated rotated geometries of the optimized geometry of one reactant/ product and translated the center of mass of the second B


Article


the Born−Oppenheimer approximation, the fact that reactants

optimized reactant/product on a sphere of radius 3 Å surrounding each of the rotated geometries of the first reactant/product. This procedure ensures uniform sampling of relative orientations of the reaction/product species. We generated 100 bimolecular orientation geometries for each set of reactants and products. Each of the reactant/product bimolecular geometries was then optimized using DFT with B3LYP, the 6-311++G** basis set, the D3 dispersion correction,34 and COSMO with dielectric constant 78.35 using the TeraChem package.22−25 The resulting geometries were then clustered using the Hungarian algorithm, which is based on the root mean square deviation between orientation geometries, and the lowest energy geometry in each cluster was selected.35 After this process, five reactants and two product geometries were selected for reaction I and four reactants and two product geometries for reaction II. In Figure 5, we show these geometries in panels (1) and (2). Given the large computational cost associated with finding MEPs, we only considered one path for reactions Iw, IIw, Ihs, IIhs, Iw+hs, and IIw+hs. The initial orientation geometries of the bimolecular reactants and products were obtained as follows. For reactions Iw, IIw, Ihs, and IIhs, a set of 10 relative orientation geometries was generated in a uniform random sampling approach, as done for reactions I and II. These geometries were subsequently optimized in TeraChem using DFT with B3LYP, the 6-311++G** basis set and the D3 correction with COSMO and dielectric constant 78.35. The lowest energy orientation geometry for each bimolecular reactant and product was chosen. For reactions Iw+hs and IIw+hs, the lowest energy orientation geometries of the reactants and products of reactions Iw, IIw, Ihs, and IIhs were selected and an additional molecule of solvent (either H2S or water) was added to their orientation geometry, which was optimized again in TeraChem. Minimum Energy Paths. An initial guess for each reaction path was obtained by interpolating between reactants and products, reactants and reactants, and products and products, using geodesic interpolation.36 The initial paths were used as input for the simplified string method21 to obtain a discrete pathway which approximates the minimum energy path and provides a guess for the location and geometry of the transition state. The Hessian of the transition state of the converged string method often had more than one minimum energy frequency and was thus a higher order saddle point. Assuming the geometry was close to the true transition state, eigenvectorfollowing37,38 was employed by using a Newton−Raphson step followed by line search, to find the first order saddle point. Indeed, by moving the transition-state geometry along all additional negative modes to reach an energy minimum in each of these directions, a first order saddle point was obtained for each path. With the new transition-state geometry, a new path was generated to connect this state to the reactant and products and the string method was applied again to get a final converged MEP with a first order saddle point transition state. In a few cases, one or two other small imaginary frequencies were still present (ω < 150 cm−1) in the eigenvalues of the Hessian of the transition state after the Newton−Raphson step, the line search, and second round of string optimization. Given their low magnitude, they were excluded from the calculation of the vibrational partition function. Microscopic Rate Constants from TST. Classical TST is based on the assumption that the reaction occurs when reactants pass through a single no return transition-state configuration.31 The approximations which lead to TST include, the validity of

are in the gas phase, at equilibrium and in the Boltzmann distribution, the assumption that there is no reflection from the barrier and no recrossing, and finally the assumption that at the transition state, which is a first order saddle point, the motion along the one-dimensional reaction coordinate can be separated from other motion and treated classically as translational motion. In general, one can show that these assumptions lead to kTST(T) ≥ kexact(T) in the purely classical limit. We take tunneling into account by using the Wigner correction, that is k(T) = κ w i g (T)·k T S T (T), where κ wig(T ) = 1 +

1 (βℏω‡)2 24

and ω⧧ = |ω⧧|, the absolute value

of the imaginary frequency of the transition state. More accurate methods and software can be employed to take tunneling into account.39−42 We account for the surrounding solvent by using the COSMO model in all electronic structure calculations (activation energies, gradients, Hessians). Some challenges in the implementation of TST include the identification of a MEP and its transition state, and the computation of accurate partition functions for the reactant and transition state. We recall that in classical TST, the expression of the rate constant, kij(T), for going from species i to species j at temperature T, in the canonical ensemble, is kij(T ) =

kBT Q‡(T ) −βEa e h Q react(T )

(1)

where β = 1/kBT with kB the Boltzmann constant, Ea = E⧧0 − Ereact 0 is the activation energy, that is the difference between the ground-state energy of the transition state, E⧧0 , and that of the ⧧ reactant ground state energy, Ereact 0 . Qreact(T), and Q (T) are the

partition functions of the reactants and transition state. Q⧧(T) is computed by excluding the reaction coordinate. The partition functions were obtained using the standard approximations (see the Supporting Information). Reaction Population Dynamics. Given N reactant species, M product species, and the (N + M)2 paths which connect these, we can write the coupled kinetic equations dC(t ;T ) = K(T )C(t ;T ) dt

(2)

where C(t;T) = {Co(t;T), ..., CNsp−1(t;T)} is the vector containing the number of molecules of each species at a given time, t and temperature, T with Nsp = N + M. We will refer to these as populations. The matrix K(T) is expressed as C


Article


Figure 3. MEPs (without the zero-point energy, ZPE, correction) for reaction I (in blue) and reaction II (in magenta). Activation-energy values include the ZPE correction. These reactions occur through two reactive events. The first starts from the lowest-energy bimolecular orientation geometry of the reactants and leads to the first product geometries (prod. 1), and the second transforms the first products into the second (prod. 2). The activation energy for reaction I is lower than that for reaction II in both steps. The reactant energy of the reaction II paths is lower than that of the reaction I paths but this is not seen as each path is rescaled by its minimum energy. The first product conformations and corresponding transition states are also indicated below the MEP plot. MEPs were computed with DFT with B3LYP, 6-311++G**, the D3 dispersion correction and COSMO with ϵ = 78.35 and no explicit solvent.

ij Nsp jj jj− ∑ k1j jj jj j ≠ 1 jj jj jj jj jj k12 jj jj jj K(T ) = jjjj jj jj k13 jj jj jj jj... jj jj jj jj jj k1N jj sp j k

k 21

... k Nsp1

k 32

...

Nsp

−

∑ k 2j j≠2

Nsp

k 23

−

∑ k 3j

...

j≠3

...

yz zz zz zz zz zz zz zz zz zz zz zz zz zz zz zz zz zz zz zz zz zz zz Nsp zz z ∑ k Nspj zzzzz z j ≠ Nsp {

k 31

...

... −

compute C(t) = U·P(t) with Pi(t) = Pi(0)·exp(λit). One can also resort to numerical integration, for example VODE28 with backward differentiation formulas. However, the K matrix is often ill conditioned, and the system of coupled ODEs is quite stiff when there is a large variation in the order of magnitude in its values. In those cases, diagonalization can lead to a large error and numerical integration packages have trouble converging. When this occurs, for a first order system, one can use exponentiation of the rate matrix, K or in general the hybrid SSA tau-leaping method described in refs30,33 to solve the coupled equations stochastically (see the Supporting Information for details of procedure). The hybrid SSA tau-leaping algorithm combines the Gillespie SSA approach,43 where at most one reaction can occur at each step, with the tau-leaping approach, where multiple reactions can occur at each time step. Therefore, it provides the advantage of a more rapid calculation of the dynamics and at the same time, a lower likelihood of unphysical negative populations during the simulation. The advantage of the hybrid SSA tau-leaping approach is evident when the coupled ODEs are not first order; indeed in those cases, the exponential is no longer the solution and an analytic solution may not be easy to find.

(3)

with kij = kij(T), the TST rate constant of the path connecting species i to j. This system of coupled first order differential equations, ODEs, is in general easy to solve exactly. One can take the exponential of the K(T) matrix to obtain C(t) = exp(K(T)·t)·C(0). Otherwise, one can diagonalize the K(T) N matrix to get its eigenvectors, U and eigenvalues, {λi}i =sp1 and D


Article


■

RESULTS AND DISCUSSION MEPs, Rate Constants, and Dynamics in the Absence of Explicit Solvent. In Figure 3, we show the MEPs obtained from the lowest-energy bimolecular reactant orientation geometry for reactions I (blue) and II (magenta) and the corresponding geometries of reactants, transitions states, and products. All MEPs for the other orientation geometries of reactants and products of reactions I and II are shown in Figure S1 of the Supporting Information. For each reactant, two initial paths were computed, one connecting to each bimolecular product orientation geometry. It was found that the second product geometry (prod. 2) is most easily obtained by going through the first (prod. 1). Therefore, the results are shown as paths connecting reactants to the first product geometries (prod. 1) and paths connecting prod. 1 to prod. 2. The ground-state energies of these reactants and products are given in Table S2 of the Supporting Information. The activation energies for reaction I are lower than those for reaction II. This is most likely due to the steric hindrance of the methyl groups in reaction II. Overall, the mechanism along each of the four paths is quite similar. As can been seen from the geometry of the transition states, the reaction mechanism is more concerted than stepwise, in fact the hydrogen atom and sulfur atoms attach to the nitrogen and carbon within two images in the first reactive event and both the H and S atoms come from the hydrosulfide anion. In the second reactive event, the NH bond rotates to transform the first product geometries into the second geometries (see prod. 1 vs prod. 2 geometries in Figure 3). In Figure 4 we show the logarithm of the rate constant as a function of 1000/T, with T, the temperature, computed by using the MEPs of the first reactive events which were shown in Figure 3. The activation energy is the main factor which influences the trend of the rate constant as a function of temperature; indeed, at

all temperatures, reaction I proceeds more rapidly than reaction II. We now consider the case when all reactant and product geometries are taken into account. In Figure 5 panels (1) and (2), we show the oriented geometries of reactants and products for reactions I and II. The population dynamics of this system (Figure 6) was obtained by solving eq 2 at T = 1000 K. We chose to compute the populations at 1000 K because the activation energy barriers are very high and thus at room temperature, the reaction rates are very low. In Figure 6, panels (1) and (2), we show the average reactant N Ci(t ), average first population, that is ⟨Creact(t )⟩ = 1/Nsp·∑i =react 1 N

product, ⟨C P1(t )⟩ = 1/Nsp·∑i =P1N

react + 1

Ci(t ), and average second N

product populations, ⟨C P2(t )⟩ = 1/Nsp·∑i =P2N + N + 1 Ci(t ). In react P1 the Supporting Information (Figure S2), we also show the populations of each reactant and product geometry, which were averaged to generate the plots in all panels of Figure 6. The initial population of each product geometry was set to zero, Ci=prod(0) = 0, and that of each reactant geometry was set to 1000, Ci=react(0) = 1000. All paths, reactant to reactant, product to product, reactant to product, and product to reactant were included and the kinetics was computed by taking the exponential of the rate matrix. In the Supporting Information, we also show that equivalent results are obtained when using integration with VODE and the hybrid SSA tau-leaping approach (Figure S2). We notice that the average populations of reactants and of products for reaction I (Figure 6, panel 1) equilibrate in a similar amount of time (Figure 6, panel 2). Sutherland et al. considered the reaction of formaldehyde with hydrogen cyanide followed by reaction I and the reaction of acetone with hydrogen cyanide followed by reaction II. They found that the second set of reactions which involve acetone is slower than the first. Experimentally,44,45 the ratio of the rate constant for the reaction of formaldehyde with cyanide with respect to that of acetone with cyanide is 8.82 × 104. Therefore, assuming a similar behavior at 1000 K, given that we do not find a large difference in our rates, we believe that the equilibration of cyanohydrin is the limiting step of the process and that reactions I and II occur on a comparable timescale. To make a quantitative comparison to future experimental results, one should also take into account association and dissociation of the reactants and eventually the role of pH.46,47 With regard to the role of including multiple pathways, we found that the rates we obtained from the half-life time of the average reactant population for reaction I and II (Figure 6) are respectively 3.5 and 33.3 times larger than those obtained from a single MEP (Table 1). This difference comes from the availability of lower reactant-to-product energy barriers (see Figure S1 in the Supporting Information). Indeed, reactants can now access these lower activation energy barriers by changing from one reactant geometry to another through a very small activation energy barrier. This emphasizes the importance of including multiple pathways in the calculation of a reaction rate constant. MEPs and Rate Constants in the Presence of Explicit Solvent. In Figure 7, we show the MEPs for reactions Iw, IIw, Ihs, IIhs, Iw+hs, and IIw+hs which include one or two molecules of the solvent explicitly. From the geometries of the transition states, shown for reaction Ihs, we see that the proton is now transferred from one

Figure 4. Logarithm of the reaction rate constant, log10 k(T), computed using TST and the Wigner correction for tunneling for reaction I (blue) and II (magenta) for the first reactive paths shown in Figure 3 as a function of 1000/T with T the temperature. The reactions with a smaller activation energy (reaction I) have a larger reaction rate constant than reaction II. As temperature decreases, this difference increases. E


Article


Figure 5. Panel (1) Reactant and product bimolecular geometries for reaction I. Panel (2) Reactant and product bimolecular geometries for reaction II. The geometries were obtained with DFT, the B3LYP functional, the 6-311++G** basis set, the D3 dispersion correction, and COSMO with dielectric constant 78.35 using the TeraChem package.

Figure 6. Panels (1,2) Average populations of reactant and product bimolecular geometries at T = 1000 K for reaction I (panel 1) and reaction II (panel 2) computed by including all paths, reactant to reactant, product to product, and reactant to product. All populations were computed by exponentiation of the rate matrix. In all panels, a rate constant is reported and was obtained by computing the half-life time from the reactant population. Given that the reactant populations are not reaching the zero-concentration limit, the half-life was obtained as the time at which the concentration is halfway between its initial value and the equilibrated non-zero final value.

Table 1. TST Rate Constants with the Wigner Correction for Tunneling,k(T), Half-Life Times, τ1/2 = ln 2/k(T), and Activation Energies,Ea, Which Include the ZPE Correction, for the MEPs of Each Reaction at 300 and 1000 Kaat 300K and 1000 K. react.

T [K]

k(T) [1/s]

τ1/2(T) [s]

Ea [kcal/mol]

react.

k(T) [1/s]

τ1/2(T) [s]

Ea [kcal/mol]

I (1) Iw Ihs Iw+hs

300

3.31 × 10−20 4.01 × 10−5 2.45 × 10−1 6.76 × 10−3

2.10 × 1019 1.73 × 104 2.83 × 100 1.03 × 102

45.17 24.40 19.50 21.71

II (1) IIw IIhs IIw+hs

1.95 × 10−23 2.70 × 10−7 1.13 × 10−3 7.40 × 10−4

3.55 × 1022 2.56 × 106 6.11 × 102 9.37 × 102

50.12 28.16 23.13 23.52

I (1) Iw Ihs Iw+hs

1000

5.81 × 103 3.24 × 108 4.60 × 109 5.94 × 109

1.19 × 10−4 2.14 × 10−9 1.51 × 10−10 1.17 × 10−10

45.17 24.40 19.50 21.71

II (1) IIw IIhs IIw+hs

5.02 × 102 4.95 × 107 4.28 × 109 4.08 × 109

1.38 × 10−3 1.40 × 10−8 1.62 × 10−10 1.70 × 10−10

50.12 28.16 23.13 23.52

The number (1) in parenthesis in the react. column indicates the first reactive event for reactions I and II.

a

of the solvent molecules and, within a couple of steps, the HS− anion attaches to form a bond with the carbon atom. This process, seen in all reactions, is energetically favorable (Table 1). Indeed, adding H2S, H2O, or both leads to a much lower activation-energy barrier, ΔEa ≈ 20 kcal/mol. Overall, when

looking at the two reactive events, we find that the two-step process is energetically more favorable when H2S is the proton donor. This is most likely because the pKa of H2S is less than that of water (pKa(H2S) ≈ 7, pKa(H2O) ≈ 14) and thus it is much more likely for H2S to dissociate and donate a proton. From the MEPs, F


Article


more acidic than water, as it has a pKa of about 7 with respect to the pKa of water which is 14. When the solvent was included explicitly (water and H2S), we also found that reaction I proceeds more rapidly than reaction II by 1 order of magnitude at 300 K, which indicates that in the reaction of formaldehyde/acetone with hydrogen cyanide (which we do not study here) followed by the reaction of formalcyanohydrin/acetone cyanohydrin with hydrogen sulfide, the first reaction is the limiting step, as it is much slower for acetone.44 With this work, we have confirmed the importance of including various reactant and product orientation geometries, as well as explicit solvent molecules to determine the kinetics accurately. Finally, we have seen that molecular reactivity in the prebiotic environment is quite different when multiple reaction paths are accounted for. Knowing this will be helpful in a future study of the full network of prebiotic reactions.

■

ASSOCIATED CONTENT

S Supporting Information *

Figure 7. MEPs for reaction Iw, Ihs, and Iw+hs (in blue) and reaction IIw, IIhs and IIw+hs (in magenta) without the zero-point energy, ZPE, correction. Activation energy values include the ZPE correction. MEPs were computed with DFT with B3LYP, 6311++G**, the D3 dispersion correction, and COSMO with ϵ = 78.35 and no explicit solvent. The energies are rescaled with respect to the minimum energy of the reactive paths. Geometries of the reactants, transitions states, and products for reaction Ihs are shown around the MEP plots; for all other pathways, the geometries can be found in the Supporting Information.

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.9b05735. Description of the calculation of the partition functions and of the hybrid SSA tau-leaping approach. Energies of the reactants and products in reactions I and IIall reactant to product MEPs for reactions I and II, geometries along the minimum energy pathways when the solvent is included explicitly, and comparison of the rate constant values when computed with TST or with TST and the Wigner correction for tunneling (PDF)

■

we computed the rate constants and half-life times (Table 1) and found that reactions Iw, Ihs, and Iw+hs have a larger rate constant than reactions IIw, IIhs, and IIw+hs as had been found without including the solvent explicitly.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (S.V.). *E-mail: [email protected] (T.J.M.).

■

CONCLUSIONS From our study of reaction I with formalcyanohydrin and reaction II with acetone cyanohydrin in the absence of an explicit solvent, we found that including multiple reactants and product geometries leads to changes in the kinetics. In fact, when all paths were included, the rate constant computed at 1000 K was about 4 to 33 times larger than the one obtained by including a single MEP from the most stable reactant orientation. When comparing reaction I to reaction II, the reaction mechanism was quite similar. The hydrogen and sulfur atoms came from the same reactant molecule and both bonded within a couple of images to the nitrogen and carbon atoms of the cyanohydrins. We also observed that reaction I proceeds more rapidly than reaction II. We then studied the role of the solvent. We saw that, for all reactions, the addition of explicit water and/or hydrogen sulfide to the reactants leads to a lowering of the activation-energy barrier of over 20 kcal/mol with respect to the case when only the implicit COSMO model is used. Furthermore, a different reaction mechanism is observed in the presence of an explicit solvent. The proton is now obtained from the solvent molecule and the sulfur from the HS− anion. When comparing the role of H2S with respect to that of H2O as a proton donor, we did not see any significant change in the mechanism. However, when hydrogen sulfide is the proton donor, the process is energetically more favorable. This can be explained by the fact that H2S is

ORCID

Stéphanie Valleau: 0000-0003-0499-2054 Todd J. Martínez: 0000-0002-4798-8947 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS S.V. would like to thank John Sutherland for providing his feedback on the computational kinetics results as well as information on his experimental work. S.V. would also like to thank Xiaolei Zhu for useful discussions on finding the MEP and transition states, David Sanchez for his code to filter reactant and product geometries using the Hungarian algorithm, and Rafał Szabla for discussions on dispersion in DFT calculations of MEPs. This work was supported the Simons Foundation Simons Collaboration on the Origins of Life (S.V., SCOL, 384805) and by the Office of Naval Research (N00014-17-1-2875).

■

REFERENCES

(1) Orgel, L. E. Prebiotic chemistry and the origin of the RNA world. Crit. Rev. Biochem. Mol. Biol. 2004, 39, 99−123. (2) Mullen, L. B.; Sutherland, J. D. Simultaneous nucleotide activation and synthesis of amino acid amides by a potentially prebiotic multicomponent reaction. Angew. Chem., Int. Ed. 2007, 46, 8063−8066. (3) Cleaves, H. J.; Chalmers, J. H.; Lazcano, A.; Miller, S. L.; Bada, J. L. A reassessment of prebiotic organic synthesis in neutral planetary atmospheres. Origins Life Evol. Biospheres 2008, 38, 105−115.

G


Article

The Journal of Physical Chemistry A (4) McCollom, T. M. Miller-Urey and beyond: what have we learned about prebiotic organic synthesis reactions in the past 60 years? Annu. Rev. Earth Planet. Sci. 2013, 41, 207−229. (5) Islam, S.; Bučar, D.-K.; Powner, M. W. Prebiotic selection and assembly of proteinogenic amino acids and natural nucleotides from complex mixtures. Nat. Chem. 2017, 9, 584−589. (6) Patel, B. H.; Percivalle, C.; Ritson, D. J.; Duffy, C. D.; Sutherland, J. D. Common origins of RNA, protein and lipid precursors in a cyanosulfidic protometabolism. Nat. Chem. 2015, 7, 301−307. (7) Klamt, A.; Schüürmann, G. COSMO: a new approach to dielectric screening in solvents with explicit expressions for the screening energy and its gradient. J. Chem. Soc., Perkin Trans. 2 1993, 799−805. (8) Cramer, C. J.; Truhlar, D. G. Implicit solvation models: equilibria, structure, spectra and dynamics. Chem. Rev. 1999, 99, 2161−2200. (9) Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. Adding explicit solvent molecules to continuum solvent calculations for the calculation of aqueous acid dissociation constants. J. Phys. Chem. A 2006, 110, 2493− 2499. (10) De Sterck, B.; Vaneerdeweg, R.; Du Prez, F.; Waroquier, M.; Van Speybroeck, V. Solvent effects on free radical polymerization reactions: the influence of water on the propagation rate of acrylamide and methacrylamide. Macromolecules 2010, 43, 827−836. (11) Steinfeld, J. I.; Francisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics; Prentice Hall Englewood Cliffs: New Jersey, 1989. (12) Peters, B. Reaction Rate Theory and Rare Events; Elsevier Science, 2017. (13) Evans, M. G.; Polanyi, M. Some applications of the transition state method to the calculation of reaction velocities, especially in solution. Trans. Faraday Soc. 1935, 31, 875−894. (14) Wigner, E. The transition state method. Trans. Faraday Soc. 1938, 34, 29−41. (15) Truhlar, D. G.; Garrett, B. C. Variational transition-state theory. Acc. Chem. Res. 1980, 13, 440−448. (16) Mandrà, S.; Valleau, S.; Ceotto, M. Deep nuclear resonant tunneling thermal rate constant calculations. Int. J. Quantum Chem. 2013, 113, 1722−1734. (17) Miller, W. H.; Zhao, Y.; Ceotto, M.; Yang, S. Quantum instanton approximation for thermal rate constants of chemical reactions. J. Chem. Phys. 2003, 119, 1329−1342. (18) Yamamoto, T.; Miller, W. H. On the efficient path integral evaluation of thermal rate constants within the quantum instanton approximation. J. Chem. Phys. 2004, 120, 3086−3099. (19) Miller, W. H.; Schwartz, S. D.; Tromp, J. W. Quantum mechanical rate constants for bimolecular reactions. J. Chem. Phys. 1983, 79, 4889−4898. (20) Wigner, E. On the quantum correction for thermodynamic equilibrium. Phys. Rev. 1932, 40, 749−759. (21) E, W.; Ren, W.; Vanden-Eijnden, E. Simplified and improved string method for computing the minimum energy paths in barriercrossing events. J. Chem. Phys. 2007, 126, 164103. (22) Titov, A. V.; Ufimtsev, I. S.; Luehr, N.; Martinez, T. J. Generating efficient quantum chemistry codes for novel architectures. J. Chem. Theory Comput. 2012, 9, 213−221. (23) Hohenstein, E. G.; Bouduban, M. E. F.; Song, C.; Luehr, N.; Ufimtsev, I. S.; Martínez, T. J. Analytic first derivatives of floating occupation molecular orbital-complete active space configuration interaction on graphical processing units. J. Chem. Phys. 2015, 143, 014111. (24) Isborn, C. M.; Luehr, N.; Ufimtsev, I. S.; Martínez, T. J. Excitedstate electronic structure with configuration interaction singles and Tamm−Dancoff time-dependent density functional theory on graphical processing units. J. Chem. Theory Comput. 2011, 7, 1814−1823. (25) Kulik, H. J.; Luehr, N.; Ufimtsev, I. S.; Martinez, T. J. Ab Initio quantum chemistry for protein structures. J. Phys. Chem. B 2012, 116, 12501−12509. (26) Jónsson, H.; Mills, G.; Jacobsen, K. W. Classical and Quantum Dynamics in Condensed Phase Simulations; World Scientific, 1998; pp 385−404.

(27) Peters, B.; Heyden, A.; Bell, A. T.; Chakraborty, A. A growing string method for determining transition states: comparison to the nudged elastic band and string methods. J. Chem. Phys. 2004, 120, 7877−7886. (28) Cohen, S.; Hindmarsh, A. CVODE, a stiff/nonstiff ODE solver. Computers in Physics; AIP Publishing, 1996; pp 138−143. (29) Gillespie, D. T. Stochastic simulation of chemical kinetics. Annu. Rev. Phys. Chem. 2007, 58, 35−55. (30) Cao, Y.; Gillespie, D. T.; Petzold, L. R. Efficient step size selection for the tau-leaping simulation method. J. Chem. Phys. 2006, 124, 044109. (31) Pechukas, P. Transition state theory. Annu. Rev. Phys. Chem. 1981, 32, 159−177. (32) Sheppard, D.; Terrell, R.; Henkelman, G. Optimization methods for finding minimum energy paths. J. Chem. Phys. 2008, 128, 134106. (33) Cao, Y.; Samuels, D. C. Discrete stochastic simulation methods for chemically reacting systems Methods in Enzymology; Elsevier, 2009; Vol. 454, pp 115−140. (34) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density function dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (35) Banerjee, S.; Liu, F.; Sanchez, D. M.; Martínez, T. J.; Zare, R. N. Pomeranz-Fritsch synthesis of isoquinoline: gas-phase collisional activation opens additional reaction pathways. J. Am. Chem. Soc. 2017, 139, 14352−14355. (36) Zhu, X.; Thompson, K. C.; Martínez, T. J. Geodesic interpolation for reaction pathways. J. Chem. Phys. 2019, 150, 164103. (37) Cerjan, C. J.; Miller, W. H. On finding transition states. J. Chem. Phys. 1981, 75, 2800−2806. (38) Wales, D. J. Rearrangements of 55-atom Lennard-Jones and (C60)55 clusters. J. Chem. Phys. 1994, 101, 3750−3762. (39) Barker, J. R.; Nguyen, T. L.; Stanton, J. F.; Aieta, C.; Ceotto, M.; Gabas, F.; Kumar, T. J. D.; Li, C. G. L.; Lohr, L. L.; Maranzana, A.; Ortiz, N. F.; Preses, J. M.; Simmie, J. M.; Sonk, J. A.; Stimac, P. J. MultiWell-2019 Software Suite; Barker J. R., University of Michigan: Ann Arbor, Michigan, USA, 2019; http://clasp-research.engin.umich. edu/multiwell/. (40) Garrett, B. C.; Truhlar, D. G. Semiclassical tunneling calculations. J. Phys. Chem. 1979, 83, 2921−2926. (41) Conte, R.; Aspuru-Guzik, A.; Ceotto, M. Reproducing deep tunneling splittings, resonances, and quantum frequencies in vibrational spectra from a handful of direct ab initio semiclassical trajectories. J. Phys. Chem. Lett. 2013, 4, 3407−3412. (42) Peters, B.; Bell, A. T.; Chakraborty, A. Rate constants from the reaction path Hamiltonian. II. Nonseparable semiclassical transition state theory. J. Chem. Phys. 2004, 121, 4461−4466. (43) Gillespie, D. T. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 1977, 81, 2340−2361. (44) Schlesinger, G.; Miller, S. L. Equilibrium and kinetics of glyconitrile formation in aqueous solution. J. Am. Chem. Soc. 1973, 95, 3729−3735. (45) Guthrie, J. P. Rate-equilibrium correlations for the aldol condensation: an analysis in terms of Marcus theory. J. Am. Chem. Soc. 1998, 120, 1688−1694. (46) Hwang, T.; Goldsmith, B. R.; Peters, B.; Scott, S. L. Watercatalyzed activation of H2O2 by methyltrioxorhenium: a combined computational−experimental study. Inorg. Chem. 2013, 52, 13904− 13917. (47) Goldsmith, B. R.; Hwang, T.; Seritan, S.; Peters, B.; Scott, S. L. Rate-enhancing roles of water molecules in methyltrioxorheniumcatalyzed olefin epoxidation by hydrogen peroxide. J. Am. Chem. Soc. 2015, 137, 9604−9616.

H


Reaction Dynamics of Cyanohydrins with ... - ACS Publications

Recommend Documents