Jonny
Proppe
,
Tamara
Husch
,
Gregor N.
Simm
and
Markus
Reiher
*
Laboratory of Physical Chemistry, ETH Zürich, Zürich, Switzerland. E-mail: markus.reiher@phys.chem.ethz.ch; Fax: +41 446331594; Tel: +41 446334308
First published on 6th July 2016
For the quantitative understanding of complex chemical reaction mechanisms, it is, in general, necessary to accurately determine the corresponding free energy surface and to solve the resulting continuous-time reaction rate equations for a continuous state space. For a general (complex) reaction network, it is computationally hard to fulfill these two requirements. However, it is possible to approximately address these challenges in a physically consistent way. On the one hand, it may be sufficient to consider approximate free energies if a reliable uncertainty measure can be provided. On the other hand, a highly resolved time evolution may not be necessary to still determine quantitative fluxes in a reaction network if one is interested in specific time scales. In this paper, we present discrete-time kinetic simulations in discrete state space taking free energy uncertainties into account. The method builds upon thermo-chemical data obtained from electronic structure calculations in a condensed-phase model. Our kinetic approach supports the analysis of general reaction networks spanning multiple time scales, which is here demonstrated for the example of the formose reaction. An important application of our approach is the detection of regions in a reaction network which require further investigation, given the uncertainties introduced by both approximate electronic structure methods and kinetic models. Such cases can then be studied in greater detail with more sophisticated first-principles calculations and kinetic simulations.
Quantum chemical methods are well suited for describing energy changes due to changes in the electronic structure of reacting molecules if these electronic effects govern the overall energetics of the process. Usually, structures considered relevant as stable intermediates or transition states are optimized and their energies are compared to identify the relevant reaction paths. Clearly, this approach is limited, especially if carried out manually, to a rather small number of structures only. For predictive work on systems for which little or no experimental information is known, the exploration of potentially important structures becomes an immense task. Several approaches exist to overcome this issue. In reactive molecular dynamics simulations,^{9–17} for example, the nuclear equations of motion are solved to explore and sample configuration space. By contrast, heuristics-guided exploration approaches are based on rules derived from chemical concepts.^{18–22} By applying predefined (possibly alchemical) transformation rules to create new chemical species, explorations in configuration space are greatly accelerated. Recently, we proposed a fully automated heuristics-guided exploration protocol^{22} in which the heuristic rules rest on reactivity descriptors derived from quantum mechanics.
It is important to understand that to theoretically grasp the kinetics of complex reaction networks, we must be prepared to investigate an enormous number of possible intermediates (on different potential energy surfaces) not generated by simple conformational changes but by the sheer number of chemically different reactants. For truly complex chemical reaction networks, no universal protocol based on quantum chemical calculations has been established so far that would span the whole range of steps from molecular and electronic structure optimization to detailed kinetic modeling. However, significant progress in all research areas that would contribute to the establishment of such a protocol has already been made. Given the algorithmic and hardware developments accomplished in the past two decades, it should be feasible to establish such a protocol in a single, integrated implementation.
Clearly, various choices and approximations need to be made and hence the protocol to be established will not be unique. Still, we demand the development of such a protocol be subjected to constraints that will make it universally applicable. Besides, we are faced with the fact that quantum chemical raw data are affected by method-inherent errors and need to be augmented by nuclear motion and temperature corrections before they can be subjected to kinetic modeling. Hence, if we must be prepared to make certain assumptions and approximations, we expect from our protocol that the violation of an approximation can be identified within the protocol and overcome by approaches beyond the realm of the protocol's standard methods. This way, we may be able to identify possible breaches that point to more sophisticated theoretical approaches to be applied. If the number of such breaches is small, then the general basis of the protocol, which will be quantum chemical methods in our case, will remain valid. And clearly, in view of the successes of quantum chemical reaction mechanism elucidation, we have good reason to believe that this is possible. Obviously, this will only be possible if we have uncertainty measures at hand that allow us to assess the accuracy of individual simulation results. For example, Vlachos and co-workers recently applied Bayesian statistics to predict rate constants for chemical processes on surfaces.^{23}
The ingredients of a general protocol for the generation and analysis of chemical reaction networks are: (1) the automated exploration of possibly relevant intermediates and transition states, (2) the determination of free energies for reactions in condensed phase, (3) systematic error estimation based on, for instance, Bayesian statistics, and (4) the kinetic modeling of the reaction network emerging.
In this work, we apply components (2)–(4) of the general protocol to a complex chemical reaction network in aqueous solution: the formose reaction. It is our goal to establish protocol-inherent validation measures that keep track of the validity of the assumptions made and that may point to advanced theoretical approaches to deliver more reliable data if needed. Moreover, our analysis is intended to be a general feasibility analysis of this protocol that will, as we shall show, point to interesting future developments.
As experimental kinetics can only examine a limited number of chemical species, thorough theoretical kinetic models corresponding to complex reaction networks spanning several time scales are desired. For the construction of a general-purpose (mass-action) kinetic model, rate constants are the essential elements to be determined.
Conventional transition state theory (TST) provides a simple approach to calculate rate constants for isothermal reactions. It is assumed in conventional TST that a reaction coordinate along a Born–Oppenheimer potential energy surface is orthogonally intersected by a hyperplane in such a way that once crossed by a trajectory starting from a reactant state, that trajectory ends in the corresponding product state.^{44} This crossing point is approximated by the first-order saddle point of a reaction coordinate. Given a canonical ensemble of microstates, for which the number of molecules N, the temperature T, and the volume V are constant, the thermal rate constant k(T) of a reaction from a reactant to a product crossing the corresponding transition state depends on the Helmholtz free energy difference between reactant and transition state, ΔA^{‡,*}, through an exponential function,^{45}
(1) |
The only quantity we need to determine for the construction of a kinetic model based on conventional TST is the Helmholtz free energy A for all intermediates and transition states contained in a given reaction network. A is determined by the, in our case canonical, partition function Q = Q(N,V,T) through A(T,Q) = −k_{B}TlnQ, where all energy states of the system of N molecules enter Q.
The direct evaluation of Q is in general not feasible. However, a well-established procedure exists^{50} to approximate gas-phase free energies. In the gas phase, Q may be factorized into a product of N molecular partition functions q of N isolated, indistinguishable, and non-interacting molecules, Q = q^{N}/N!. This factorization enables the calculation of A based on an isolated-molecule quantum chemical calculation. A results from the internal energy U and the temperature-weighted entropy S,
A(T,Q) = U(T,Q) − TS(T,Q). | (2) |
For our set-up and standard state, U can be decomposed into the sum of the temperature-independent electronic energy, E_{elec}, the zero-point vibrational energy, ZPE, and remaining thermal contributions to the internal energy U_{rest}(T,Q),
A(T,Q) = N_{A}E_{elec} + N_{A}ZPE + U_{rest}(T,Q) − TS(T,Q), | (3) |
The molecular partition function q is usually approximated by a product of electronic (elec), q_{elec}, translational (trans), q_{trans}, rotational (rot), q_{rot}, and vibrational (vib), q_{vib}, contributions – deliberately neglecting the coupling of the degrees of freedom (such as rovibrational coupling).
For the evaluation of q_{elec}, we may assume that electronically excited states are high in energy and cannot be excited thermally at a given temperature.^{50} However, for intermediates with small HOMO–LUMO gaps, this may be taken as an indication for a necessary extension of this standard approach. Hence, only spin and orbital or point-group degeneracy needs to be taken into account.
To evaluate q_{trans}, the particle in a box model is employed to determine the energy states associated with the translation of the center of mass of the molecule in V.^{50}q_{trans} depends on the mass of the molecule and V. If we neglect the existence of isotopomers, the mass is easily calculated for the most abundant isotopes. Recall that V is the volume for the chosen standard state.
To evaluate q_{rot}, the molecule is treated as a rigid rotor.^{50} The molecular principal moments of inertia and the symmetry number of the molecule enter q_{rot}. The error introduced by this assumption is severe, for instance, for microsolvated molecules, which feature not only internal rotational degrees of freedom, but also their coupling with the external rotations, i.e., the moments of inertia depend on internal rotational degrees of freedom.^{51}
q _{vib} is usually approximated by all energy eigenvalues of harmonic oscillators that describe the 3m − 6 normal coordinates of an m-atomic non-linear molecule.^{50} It is well known that the harmonic approximation will break down, if the potential energy surface deviates significantly from the quadratic harmonic potential,^{50,52–56} for instance, for highly anharmonic modes, which are better described as internal rotations.^{57} The evaluation of the contribution of these anharmonic frequencies is computationally expensive,^{58–63} and therefore, not feasible for the exploration of a large reaction network. If, however, reactants and products exhibit similar harmonic vibrational frequencies, error cancellation will occur (see, e.g., ref. 64).
The standard-state free reaction energy in solution, , for a reaction R → P in solution can be obtained as the sum of the standard-state free reaction energy in the gas phase, , and the difference of the free energies of solvation of the reagents, , when employing a thermodynamic cycle as illustrated in Fig. 1.
Fig. 1 Thermodynamic cycle for the calculation of the standard-state Helmholtz free reaction energy for the reaction R → P in solution (‘sol’) from gas-phase (‘gas’) and solvation (‘solv’) data. |
Implicit solvation models have been specifically devised to calculate , without modeling the solvent explicitly.^{65–67} The solute is immersed into a cavity in the continuum leading to parameter-dependent solvent–solute interactions. Implicit solvent models vary in their description of these interactions between the solvent and the solute.^{66,68–70} The parameters in the implicit solvation models are determined through parametrization to experimental Gibbs free energies of solvation.^{66} The Gibbs free energy of solvation is taken to be equal to the Helmholtz free energy of solvation,^{66} which neglects volume changes. This effect was shown to be small for several test cases.^{71} Here, we assume that all effects associated with the transfer of the molecule from an ideal gas phase to the solution phase are absorbed in .^{57,72} This assumption is reasonable for small, rigid molecules, whose structures do not significantly change upon solvation. For such molecules, modern implicit solvation models (e.g., SM12)^{70} exhibit a mean unsigned error of about 2 kJ mol^{−1} in comparison to experimental data.^{70,72} A more pronounced error can be expected for charged species.^{70,72,73} The error can, however, be reduced when adding explicit solvent molecules and averaging conformationally.^{74–76}
The choice of implicit solvent models limits the choice of possible temperatures, because they are usually parametrized to Gibbs free energies of solvation at ambient temperature.^{66} The fit to experimental data at one temperature implies that a subdivision of the into and is not possible. Note, however, that Chamberlin et al. introduced a temperature-dependent implicit solvation model, which could expand the range of accessible temperatures.^{77}
Furthermore, many popular continuum solvation models assume that thermal equilibrium between solute and solvent is reached instantaneously. This may be inadequate for reactive intermediates.^{67,78}
In this work, a reaction or activation free energy in solution, , is obtained from the difference of the electronic energies, ΔE_{elec}, the difference of the zero-point vibrational energies, ΔZPE, the difference of the thermal contributions to the free energy in the gas phase, , and the difference of the free energies of solvation, ,
(4) |
Here, Q_{solv} is a place holder for the continuum modeling of . We highlight the separation of contributions that do not depend on a partition function (ΔE_{elec}, ΔZPE) from those that do depend on it (, ). Recently, we have demonstrated how the error associated with ΔE_{elec} can be assessed by applying Bayesian statistics^{79} (see also Section 4). This error may be considered a lower bound for the error of if the error of the other contributions is neglected. For the error estimation of , ZPE, and the individual contributions must be investigated.
The model-inherent errors in Q, which is employed to calculate , are difficult to evaluate. We may assume that the error of q_{vib} is dominant. As the harmonic approximation is a severe approximation for low-frequency modes,^{52–56} the anharmonic q_{vib} is required to assess the effect of this approximation. The application of scaling factors for harmonic frequencies^{62,80} is, however, not sufficient to obtain a reference anharmonic q_{vib}, because scaling factors neither correct the form of the quadratic harmonic potential nor the equidistance of the energy levels. Recently, procedures were outlined how the full-dimensional potential energy surface can be dissected as a sum of independent one-dimensional potentials for each vibrational mode.^{53–56} The one-dimensional potentials are sampled along the normal coordinates. The energy levels of the system are obtained by solving the one-dimensional nuclear Schrödinger equations for each mode.^{53–56} Hence, the deviation of the harmonic ZPE from the anharmonic ZPE is readily obtained. It is then possible to assess the error of the harmonic q_{vib} by comparison with the anharmonic q_{vib} obtained by explicit summation over all vibrational modes. This error can be considered a lower bound for the error of since mode–mode coupling effects and errors due to the approximations in q_{trans} and q_{rot} are not taken into account.
Accurate, theoretical reference data are difficult to obtain for .^{66,81} Alternatively, the error of can be estimated by comparison to available experimental data for benchmark sets (e.g., the Minnesota Solvation Database^{82}), assuming transferability.
In general, it is difficult to predict the error of density functional calculations.^{83} To overcome this issue, Jacobsen, Sethna, Nørskov, and co-workers devised a scheme for systematic error estimation of DFT results based on non-hybrid density functionals.^{85–88} By generating an ensemble of exchange–correlation functionals, a mean and a variance could be assigned to each result (see also ref. 89 and 90).
Based on the work of Jacobsen, Sethna, Nørskov, we developed a novel approach for the construction of reliable density functionals with Bayesian error estimation capabilities.^{79} Our ansatz was tailored to overcome the transferability problem by relying on system-focused reference data and to exploit the better accuracy of hybrid functionals.
Instead of considering only the best-fit parameter (as commonly done with standard exchange–correlation functionals), we assign a conditional probability distribution to a linear parameter a in the exchange–correlation functional,
(5) |
(6) |
(7) |
By a system-focused reparametrization of LC-PBE0, the long-range corrected (LC) version of the density functional PBE0,^{91–93} we were able to reliably estimate errors of calculated reaction energies.^{79} Hereinafter, we refer to such a functional as LC*-PBE0 in order to emphasize that the original parameters were modified (in this work according to data related to the formose reaction). In our previous study,^{79} we concluded that four parameters in this exchange–correlation functional need to be modified to achieve accurate relative energies and reliable error estimates for a specific chemical system.
For an accurate reparametrization, the reference dataset needs to be representative for the system to be studied. In this study, contains structures of intermediates and transition states of the formose reaction. Specifically, = {(A_{i},B_{i})} consists of pairs of structures on the same potential energy surface, i.e., structures with the same number and type of atomic nuclei, the same number of electrons, and the same electronic spin state. Then, the electronic energy difference ΔE_{elec,i} = E_{elec}(B_{i}) − E_{elec}(A_{i}) between the two structures A_{i} and B_{i} of data set entry i defines the cost function C,
(8) |
To assess the transferability of the reparametrized functional, the dataset was arbitrarily split into a training set and a test set with 25 and 17 entries, respectively. By minimizing C with respect to the training set with the L-BFGS-B algorithm,^{94} a new set of parameter values for the LC*-PBE0 functional was obtained (see ESI†).
In Tables 1 and 2, the accuracy of LC*-PBE0 (in comparison to standard functionals) with respect to the training set and test set is given. It can be seen that LC*-PBE0 is significantly more accurate than most standard functionals considered here. The optimized parameters of LC*-PBE0 are close to those of PBE0 (see ESI†), which explains why the functionals are of similar accuracy. Due to its additional parameters, and therefore, higher flexibility LC-PBE0 was chosen over PBE0 for the re-parametrization. Nonetheless, with a largest absolute deviation between 8–10 kJ mol^{−1}, it is clear that error estimation is still necessary.
LAD | MAD | MSD | |
---|---|---|---|
B3LYP | 18.7 | 7.6 | 1.7 |
B3LYP-D3 | 22.2 | 7.0 | 1.6 |
BP86 | 28.5 | 7.3 | 1.8 |
BP86-D3 | 32.6 | 6.5 | 1.7 |
LC-PBE0 | 37.2 | 13.6 | 0.7 |
M06-2X | 20.9 | 7.5 | 1.2 |
M06-2X-D3 | 20.8 | 7.4 | 1.2 |
M06-L | 19.4 | 9.6 | 2.1 |
M06-L-D3 | 19.5 | 9.6 | 2.1 |
PBE | 28.8 | 6.2 | 1.6 |
PBE0 | 13.5 | 5.7 | 1.2 |
PBE0-D3 | 16.3 | 5.2 | 1.2 |
TPSS | 37.3 | 14.3 | 3.4 |
TPSS-D3 | 33.2 | 13.9 | 3.3 |
TPSSh | 32.3 | 13.6 | 3.0 |
TPSSh-D3 | 29.2 | 13.1 | 2.9 |
LC*-PBE0 | 9.8 | 3.7 | 1.0 |
LAD | MAD | MSD | |
---|---|---|---|
B3LYP | 14.7 | 6.0 | −0.1 |
B3LYP-D3 | 20.0 | 6.4 | 0.8 |
BP86 | 19.6 | 6.7 | 0.4 |
BP86-D3 | 25.0 | 7.8 | 1.5 |
LC-PBE0 | 27.5 | 8.4 | −1.1 |
M06-2X | 12.0 | 4.7 | −0.1 |
M06-2X-D3 | 12.0 | 4.7 | −0.1 |
M06-L | 20.0 | 7.5 | 1.5 |
M06-L-D3 | 20.3 | 7.6 | 1.5 |
PBE | 19.9 | 6.4 | 0.7 |
PBE0 | 11.9 | 4.0 | −0.1 |
PBE0-D3 | 14.7 | 4.0 | 0.5 |
TPSS | 16.6 | 6.4 | −1.0 |
TPSS-D3 | 17.6 | 7.4 | −0.3 |
TPSSh | 15.0 | 5.4 | −1.2 |
TPSSh-D3 | 15.4 | 6.2 | −0.4 |
LC*-PBE0 | 8.0 | 2.7 | 0.1 |
In Fig. 2 and 3, LC*-PBE0 (with ±σ error bars, calculated from an ensemble of N_{BEE} = 50 functionals as described in the ESI†) is compared to contemporary density functionals with respect to the training set and test set, respectively.
Fig. 2 Errors of LC*-PBE0 (with error bars indicating ±σ) and several standard functionals with respect to the training set (D1–D25). The dashed lines indicate an error of ±4.2 kJ mol^{−1}. |
Fig. 3 Errors of LC*-PBE0 (with error bars indicating ±σ) and several standard functionals with respect to the test set (D26–D42). The dashed lines indicate an error of ±4.2 kJ mol^{−1}. |
For both data sets, we observe that the error is at least within ±4.2 kJ mol^{−1} (≈1 kcal mol^{−1}), unless the error estimate reported by the functional indicates otherwise (i.e., σ > 4.2 kJ mol^{−1}). While there are relative energies for which the errors are underestimated (D2, D4, and D25 in the training set and D30 and D38 in the test set), considering the diversity of this reference set and the error of some standard functionals (see also Tables 1 and 2), the accuracy of the error estimation is satisfactory.
Since differential equations describing chemical processes are generally coupled, analytical integration becomes rapidly impossible. Therefore, numerical integration is the standard method of choice. However, given that reaction networks may be arbitrarily large and entangled, numerical integration may become inefficient,^{95} especially if the underlying process spans multiple time scales. For this purpose, a variety of approaches were designed focusing on the simplification of kinetic models.^{96}
We assign time-dependent population densities y_{i}(t) with i ∈ {1,…, n} (here: molar concentrations) to the vertices, and rate constants k^{forward}_{j} and k^{backward}_{j} with j ∈ {1,…, l} to the edges. Assuming detailed balance, the rate of a reaction pair f_{j} reads
(9) |
Note that we only consider chemical reactions with a molecularity smaller than or equal to 2. We consider this to be a good assumption for solution chemistry as long as solvent–solvent reactions are unlikely to occur (e.g., for the formose reaction in pure formaldehyde, a trimolecular reaction of formaldehyde to 1,3,5-trioxane could be feasible^{37}).
(10) |
Our objective is to integrate this kinetic model such that concentration trajectories are obtained from the initial conditions (feed of reactants, temperature) up to thermodynamic equilibrium.
For the development of our kinetic simulation algorithm, we were inspired by two such simplification approaches, namely Markov State Models (MSMs)^{97,98} and Computational Singular Perturbation (CSP).^{99,100} MSMs were developed for molecular dynamics simulations, where the phase space is decomposed into microstates such that a formerly continuous trajectory becomes a jump process, which is no longer Markovian (memoryless). Since local information is lost in a discrete phase space, the decomposition is chosen such that transitions within a microstate are much more likely to occur than transitions between microstates. This way, rapid convergence to local equilibrium can be assumed for these microstates, which recovers Markovianity. As a consequence, a kinetic model can be constructed from these discrete microstates by counting transitions between them. The microstates may in turn form macrostates (kinetic clusters) for which transitions are much more likely to occur than transitions between them. These clusters can be determined by studying the eigenvalues λ_{i} of the n × n rate matrix K = {K_{ij}}.^{101} Its elements K_{ij} are a measure of the rate for a transition from the j-th to the i-th microstate. In the case of linear kinetic models (first-order reactions only) as studied in MSMs, the rate matrix K is time-invariant and equals the Jacobian J = {J_{ij}},
(11) |
The time scale τ_{i} corresponding to the process described by the i-th eigenvalue is inversely related to the modulus of that eigenvalue,
τ_{i} = |λ_{i}|^{−1}, | (12) |
g = g_{fast} + g_{slow}. | (13) |
With this decomposition at hand, it is possible to dissipate the fast processes applying a small time step τ_{fast} in the numerical integration until g_{fast} ≈ 0. Subsequently, the slow processes can be modeled from the updated initial conditions applying a much larger time step τ_{slow}. Clearly, the larger the demanded spectral gap, the smaller the error introduced by assuming decoupling of fast and slow processes.
Since the Jacobian is time-invariant in the case of linear kinetic models, the time-scale separation is also invariant in the course of the global reaction process and needs to be examined only once. However, in non-linear kinetic models (as studied here), the Jacobian is a function of time due to the inclusion of concentrations of reaction partners.^{96} This poses a challenge to the time-scale separation as now a steady examination of the time gap is necessary to ensure valid decoupling of fast and slow processes.
One of the most robust approaches in this respect is CSP.^{102} The basis of CSP is the assumption that the concentration trajectory of a chemical process is rapidly attracted onto a slow invariant manifold Ω,^{99} which is an (n − m)-dimensional hypersurface in concentration space, where n denotes the number of species and m denotes the number of fast time scales. Consequently, τ_{m} and τ_{m+1} are the time scales of the slowest of the m fast processes and of the fastest of the (n − m) slow processes, respectively. Two subspaces, the m-dimensional subspace of fast processes and the (n − m)-dimensional subspace of slow processes, are introduced, which are spanned by m n-dimensional (column) basis vectors a_{i} (j ∈ 1,…, m) and (n − m) n-dimensional (column) basis vectors a_{j} (j ∈ m + 1,⋯, n), respectively. Furthermore, a set of n-dimensional dual (row) basis vectors b^{p} (p ∈1,⋯, n) is employed, which fulfill the condition b^{p}a_{q} = δ_{pq}, where δ_{pq} is the Kronecker delta. The decomposition ansatz for the rate vector reads
g_{fast} = [a_{1},…,a_{m}][b^{1},…,b^{m}]g, | (14) |
g_{slow} = [a_{m+1},…,a_{n}][b^{m+1},…,b^{n}]g. | (15) |
CSP approximates the basis vectors b^{p} and a_{q} by an iterative refinement procedure, where each refinement introduces an accuracy increase.^{100} The refinement procedure requires the time derivatives of the basis vectors.^{99} Therefore, computational savings due to the time-scale separation may be lost by iteratively determining the basis vectors after each time step.^{96} However, the first refinement does not involve time derivatives and already guarantees numerical stability of the simplified model.^{100}
Typically, the left and right eigenvectors of the Jacobian are chosen as an approximation for the basis vectors b^{p} and a_{q}, respectively, according to the CSP formalism. This approximation corresponds to the first refinement of the CSP basis vectors.^{95} It follows that the eigenvalues of the Jacobian can be obtained from our choice of CSP basis vectors,
λ_{i} = b^{i}Ja_{i}. | (16) |
Here, we follow an alternative approach to determine which one of the l reaction pairs contributes to the fast processes. We consider the largest modulus of eigenvalues of each of the l sub-Jacobians corresponding to the isolated reaction pairs. If a dominant modulus is larger than |λ_{m+1}|/ε (λ_{m+1} is associated with the Jacobian of the entire reaction system), the reaction pair connected to it will contribute to the fast processes. The idea behind this approach is that an eigenvalue corresponding to the entire kinetic model is approximately the sum of the eigenvalues of sub-Jacobians with similar or smaller moduli.^{95} With this approach, we introduce the assumption that a reaction pair is either included in or excluded from the fast processes, which certainly is a simplification that requires careful investigation.
Next, we propagate the fast sub-network (i.e., the sub-network containing only those edges corresponding to fast reaction pairs) to local equilibrium. The stationary distribution can be determined through a non-linear optimization algorithm^{103} or analytically for simpler networks. Due to the time-scale separation, it is assumed that this process occurs immediately, i.e., it is not resolved in the course of the kinetic simulation.
Then, the actual simulation starts. The partially equilibrated concentration vector y_{peq}(t) is propagated according to the time scale τ_{1,slow}, which corresponds to the fastest process of the Jacobian of the slow sub-network (i.e., the sub-network containing only those edges corresponding to slow reaction pairs). The update of the concentration vector reads
y(t + τ_{1,slow}) = y_{peq}(t) + g(t)τ_{1,slow}. | (17) |
After that, the Jacobian of the entire network is decomposed again to determine the fast and slow processes for the next time step.
Our kinetic simulation algorithm can be summarized as follows:
(1) Determine the number of fast time scales m by spectral decomposition of the Jacobian corresponding to the kinetic model under consideration.
(2) Identify a reaction pair as a fast one if the largest modulus of eigenvalues of its sub-Jacobian is larger than |λ_{m+1}|/ε.
(3) Propagate the fast sub-network to local equilibrium, y(t) → y_{peq}(t).
(4) Determine the time step τ_{1,slow} by decomposing the Jacobian of the slow sub-network.
(5) Update the partially equilibrated concentration vector according to eqn (17), y_{peq}(t) → y(t + τ_{1,slow}).
(6) If global equilibrium is not yet reached, repeat steps 1 to 5; otherwise, stop the simulation.
This sub-network, which already features many conceptual challenges of the entire formose reaction, is shown in Fig. 4. It represents a possible mechanism for the first steps of the formose reaction as described by Kua et al.^{37} and comprises six chemical species and five reaction pairs (ten elementary reactions Ri). We obtained all free energies in single-point calculations as described in the ESI.† In water, formaldehyde (1) is in equilibrium with its hydrated form, methanediol (2). 1 dimerizes to glycolaldehyde (3), which is a reaction with a high free activation energy (cf.Table 3). The exact mechanism of the dimerization has not been unravelled yet.^{31–34} From experimental studies it is, however, well known that the dimerization proceeds very slowly. 3 can react with water to 1,1,2-ethanetriol (5). Another possible reaction of 3 is the enolization to 1,2-ethenediol (4). The addition of 1 to 4 yields glyceraldehyde (6). This bimolecular reaction introduces a significant entanglement in the model network. The model network does not capture the autocatalytic nature of the formose reaction, in which 3 can be regenerated autocatalytically from intermediates generated from 6 in subsequent reactions.
Fig. 4 (a) Possible mechanism of the first steps in the formose reaction, (b) abstract graph representation of this reaction sub-network. |
Reactant(s) | Product(s) | ΔA^{‡,*} | k | ||
---|---|---|---|---|---|
R1 | 1 | 2 | 95.4 | 4.8 | 6.7 × 10^{−3} |
R2 | 2 | 1 | 124.9 | 13.2 | 8.1 × 10^{−10} |
R3 | 1 + 1 | 3 | 215.4 | 14.2 | 1.2 × 10^{−25} |
R4 | 3 | 1 + 1 | 311.1 | 23.0 | 1.9 × 10^{−42} |
R5 | 3 | 4 | 157.3 | 11.6 | 1.7 × 10^{−15} |
R6 | 4 | 3 | 130.8 | 10.2 | 7.5 × 10^{−11} |
R7 | 3 | 5 | 100.3 | 3.2 | 9.2 × 10^{−4} |
R8 | 5 | 3 | 119.2 | 12.3 | 8.0 × 10^{−9} |
R9 | 1 + 4 | 6 | 112.5 | 13.4 | 1.2 × 10^{−7} |
R10 | 6 | 1 + 4 | 185.4 | 23.1 | 2.0 × 10^{−20} |
Table 3 presents (standard-state Helmholtz) free activation energies (in solution), ΔA^{‡,*}, calculated according to eqn (4) and the resulting rate constants k (together with error estimates) for the reactions in the model network.
It can be seen that ΔA^{‡,*} is high (above 100 kJ mol^{−1}) for most reactions, and consequently, the reaction rates are small. In addition, most reactions have estimated errors of above 10 kJ mol^{−1}, which reflects the large uncertainty of the respective reaction rates. In Section 4, we showed that the LC*-PBE0 functional provides reliable error estimates above 4.2 kJ mol^{−1}. The estimated error for reaction R7 is below that, and therefore, most likely too small.
For the simulation we selected an absolute temperature of 298.15 K, a 1 M solution of formaldehyde in water as initial feed, and a time-gap criterion of ε = 10^{−3}. For technical details of the kinetic modeling employed here, see the ESI.†
For every set of free activation energies, it was found that all reaction pairs but (R3, R4), the dimerization of formaldehyde (1) to glycolaldehyde (3), contribute to the fast processes. Therefore, only reaction pair (R3, R4) constitutes the slow sub-network (Fig. 5).
Fig. 5 Fast (bottom left) and slow (top right) sub-networks of the reaction network shown in Fig. 4. |
The concentration trajectories from the kinetic simulation of the reaction network are shown in Fig. 6. The red curve corresponds to the trajectory obtained from the free activation energies listed in Table 3. The black curves correspond to the trajectories obtained from the free activation energies calculated from the ensemble of density functionals generated by our error estimation scheme according to eqn (6).
Fig. 6 Concentration trajectories with respect to time for chemical species 1–6 according to the reaction network shown in Fig. 4. The trajectories resulting from the free activation energies listed in Table 3 are shown in red. The other trajectories (black) result from the free activation energies calculated from the ensemble of density functionals generated by the error estimation scheme. Note that the time scale of the equilibration process is extremely large, which originates from neglecting relevant intermediates and elementary reactions in our model network. For readability reasons, all plots start after the first global time step τ_{1,slow}, which depends on the sampled free activation energies, and therefore, the onset of the trajectories is different. |
The simulated time scale of the global process exceeds the age of the universe in each case. This finding should not be interpreted in absolute terms, but it indicates that the uncatalyzed thermal formose reaction is very unlikely to occur if one starts from formaldehyde (1) alone, provided that the free activation energies and their estimated errors are reliable. It should be noted that glycolaldehyde (3) is autocatalytically regenerated in the formose reaction, which is not considered in our model network. This way, the dimerization of 1 to 3 (reaction R3) can be circumvented, which leads to an acceleration of the overall process not depending on the extremely slow reaction R3.
The concentration trajectories show clearly how sensitive rate constants are to variations in the free activation energies. For instance, the variation in time of the concentration trajectories of methanediol (2) spans almost 23 orders of magnitude (a factor of 8.7 × 10^{22} at an arbitrarily chosen concentration of y_{2} = 0.01 mol L^{−1}). Since only reactions R3 and R4 contribute to the time resolution of the chemical process, uncertainties in the corresponding free activation energies need to be responsible for this significant variation. In Table 4, properties of the fastest and slowest concentration trajectories (Fig. 6, species 2, left-most and right-most curves) are compared. For reaction R3, the free activation energy spans a range of about 60 kJ mol^{−1}, and for reaction R4, this range is about 100 kJ mol^{−1}, leading to a deviation in rate constants of about 10 and 17 orders of magnitude, respectively. Taking the concentrations of 1 and 3 (the constituents of reactions R3 and R4) at our arbitrarily chosen concentration of y_{2} = 0.01 mol L^{−1} into account, the rates of both reactions can be calculated. For both the fastest and slowest trajectories, reaction R3 is much faster than reaction R4. Therefore, we assume only reaction R3 to be relevant for the kinetic simulation. The reaction time can be roughly estimated by the inverse of the current reaction rate. In our case, the reaction time of the slowest trajectory is higher than that of the fastest trajectory by a factor of 1.4 × 10^{23}, which is quite close to the factor of 8.7 × 10^{22} determined from the concentration data of 2. Obviously, an error of this magnitude with respect to the free activation energy is far too large to quantify concentration trajectories in terms of absolute time. Moreover, it should be noted that the error introduced by choosing conventional TST to calculate rates is not considered here.
y _{1} | y _{3} | r _{R3} = k_{R3}(y_{1})^{2} | r _{R4} = k_{R4}y_{3} | |
---|---|---|---|---|
Fastest | 1.4 × 10^{−4} | 2.9 × 10^{−2} | 9.3 × 10^{−28} | 6.9 × 10^{−35} |
Slowest | 7.4 × 10^{−11} | 1.5 × 10^{−7} | 6.5 × 10^{−51} | 2.5 × 10^{−57} |
Even though the uncertainty in free activation energies strongly affects absolute time, it does not affect the qualitative flux of concentrations through the network in terms of non-crossing trajectories (Fig. 7). This finding can be explained by the distinct separation of the magnitude of the free activation energies. In Table 4, it can be seen that the free activation energy for reaction R3 in the slowest case is even lower than that for reaction R4 in the fastest case. Furthermore, the free activation energies and their uncertainties listed in Table 3 show that all reaction barriers are well separated from each other, which does not allow for an alternative reaction mechanism. Clearly, for small activation energy differences, such as found in enantioselective organocatalysis, large uncertainties would also lead to qualitatively different results.
Fig. 7 Concentration trajectories with respect to reaction progress for chemical species 1–6 according to the reaction network shown in Fig. 4. The trajectories resulting from the free activation energies listed in Table 3 are shown in red. The other trajectories (black) result from the free activation energies calculated from the ensemble of density functionals generated by the error estimation scheme. Contrary to Fig. 6, here, the trajectories are laid on top of each other. |
Qualitative validity of the kinetic simulation is also underlined by the fact that in all cases, 1,1,2-ethanetriol (5) is the main product at chemical equilibrium. The population dominance of 5 over 3 was also found experimentally by Kua et al.^{105} However, their calculated Gibbs free activation energies for the corresponding reaction pair (R7, R8)^{37} () are very similar to each other. Their Gibbs free activation energies can be directly compared to our Helmholtz free activation energies, because volume changes are neglected. Our free activation energies for the reaction pair (R7, R8) differ significantly from each other on average (). A reason for the observed difference is the choice of computational methods for the calculations (e.g., different density functional and solvation model). It might seem surprising that 5 is the main product in our simulation even though glyceraldehyde (6) is a thermodynamic sink. However, one should keep in mind that the concentration trajectory of 6 is temporarily significantly populated. To understand this finding, we need to discriminate between fast and slow processes (Fig. 5).
Considering the fast sub-network (Fig. 5), we understand that there are two unconnected channels to form 6, i.e., (1, 2) and (3, 4, 5). This picture is equivalent to reaction A + B ⇌ C, where the initial concentration difference between A and B is conserved over the course of the reaction. It follows that
Δ ≡ (y_{1} + y_{2}) − (y_{3} + y_{4} + y_{5}) | (18) |
Another feature of the conservation of Δ is that the kinetic model can be reduced to a single differential equation (see the ESI† for more details). This differential equation can be easily integrated by any conventional numeric solver. Here, we chose the standard fourth-order Runge–Kutta algorithm. We compared the result to that of our CSP-type method, where we employ an explicit Euler algorithm according to eqn (17), which is the simplest ansatz for numerical integration and known to be unstable due to the lack of an inherent time step selection. However, our CSP-type method provides the time step for the explicit Euler algorithm by continuously analyzing the Jacobian. We emphasize that both approaches to model the kinetics of the network (CSP/Euler vs. Runge–Kutta) yield identical results.
Then, the heuristics-guided exploration established for protonation reactions in a previous study^{22} needs to be extended to enable the exploration of chemical systems such as the formose reaction,^{104} which will facilitate a fully automated exploration and analysis of the formose reaction.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6fd00144k |
This journal is © The Royal Society of Chemistry 2016 |