Mean first-passage times for solvated LiCN isomerization at intermediate to high temperatures

The mean first-passage time (MFPT) is a useful estimate of an inverse reaction rate constant (for simplicity called rate in this paper) whether it is determined numerically or analytically.^1–7 For example, the inverse MFPT can be used to numerically obtain the Kramers rate³ through the direct observation of an ensemble of trajectories of a given system. Its use in resolving rates in chemical reactions has been demonstrated in several systems, including the isomerization reaction of lithium cyanide (LiNC/LiCN). The latter has received much attention,^8–20 in part, because of the early introduction of a useful model potential by Essers et al.²¹ Recently, for example, some of us²⁰ obtained the decay rate of trajectories near the normally hyperbolic invariant manifold (NHIM) for this reaction so as to demonstrate the use of nonrecrossing dividing surfaces (DSs) in time-dependent transition state theory (TST).²² 

In turn, the objective of this paper is to demonstrate a stronger connection between rates obtained using MFPT to those from flux-over-population formulas. Pollak, Grabert, and Hänggi (PGH)²³ employed the latter in constructing the theory describing the Kramers turnover²⁴ in the rates from low to high friction. The turnover was observed in the LiNC reaction²⁵ and has subsequently been seen to correlate well with the Pollak–Grabert–Hänggi (PGH) theory.^16,17,19

Here, we obtain the MFPT and the PGH theory explicitly for the backward reaction LiCN → LiNC. Perhaps surprisingly, we found that the MFPT can be useful in describing reactions of particles undergoing Brownian motion at much higher temperature than the barrier height. Throughout this work, we use the terms low, intermediate, and high to qualify our temperature range in reference to the barrier height. Intermediate temperatures are on the order but lower than the barrier height, whereas low temperatures are much lower than the barrier height by at least a factor of 10 such that the system is mostly bound and the reactions must be activated. Notably, we benchmarked the rates obtained from the MFPT to several alternative methods so as to confirm the relative accuracies across a broad range of friction and temperature.

At high temperatures, the rate of escape and the corresponding rate constant have to be carefully defined. Here, a significant proportion of the initial conditions are at initial energies above the barrier. These are not trapped and simply escape across the barrier without return (as long as an absorbing boundary condition is used on the other side). These trajectories do not contribute to a steady-state flux or the corresponding rate constant. However, there is still a small population—perhaps fleetingly so—whose energies are initially below the barrier and, hence, initially trapped. They, of course, gain energy quickly through interactions with the bath nearly instantaneously redistributing themselves into yet another Boltzmann distribution. Their rate of escape is the one that we calculate and use to determine a rate constant.

The structure of this paper is as follows. The model for the backward reaction LiCN → LiNC is summarized in Sec. II and connected to a solvent through the Langevin equation (LE) in Sec. II B. The use of PGH theory and MFPT to obtain reaction rates for this reaction is presented in Secs. II C and II D, respectively. The MFPT rates are presented in Sec. III A, and the nature of their turnover is addressed in Sec. III B. A comparison of these methods to all-atom molecular dynamics (AAMD) for the prediction of the maximum in the rates in the challenging case of high temperatures is presented in Sec. III C.

The isomerization reaction LiNC ⇌ LiCN involves the breaking and making of weak bonds between a lithium cation and a cyanide anion. LiNC is the stable conformation of the isomerization reaction, and hence, the backward reaction

is exothermic. Figure 1 shows the LiCN molecule in a Cartesian body-fixed frame.⁸ Through the approximation of separating the translation and rotation of the center of mass of the whole molecule, it is possible to reduce the equations of motion to two internal variables. These variables, R and ϑ, are known as the Jacobi coordinates: R = |R| is the distance between the NC center of mass and the lithium atom, and $ϑ=∡(rNC,R)$ is the angle between R and r_NC, where |r_NC| = r_e = 2.186a₀, and a₀ is the Bohr radius. At ϑ = 0, the potential is at a local minimum that corresponds to the metastable conformation of the LiCN reactant in the backward reaction. The global minimum of the potential lies at ϑ = π and corresponds to the stable conformation of the LiNC product.

The Hamiltonian can now be specified using the Born–Oppenheimer potential energy surface (PES)²¹ and the usual approximations^9,13 as

containing the kinetic and potential²¹ V(R, ϑ) energy in R and ϑ, where the reduced masses are $μ1=[1/mLi+1/(mC+mN)]−1$⁠, $μ2=(1/mC+1/mN)−1$⁠, and

We assume pure ⁷Li,¹²C, and ¹⁴N isotopes for the atomic masses m_Li, m_C, and m_N, respectively. The potential energy surface (PES) in Jacobi coordinates is shown in Fig. 2. The reaction saddle is marked with a cyan cross and lies at R = 4.2197 a.u. and ϑ = 0.2922π.

In this paper, we use atomic units (a.u.) for all our figures and calculated values. For example, the distances are given in terms of the Bohr radius a₀; the atomic time unit is given by ℏ/E_h, where E_h is the Hartree energy; and the mass is given in electron masses m_e. In these units, we found high-temperature rates k on the order of 10⁻⁴ a.u. (cf. Fig. 5) corresponding roughly to k = 4 × 10¹² s⁻¹ in SI units. The inverse of such a rate is comparable to the time scale of the isomerization dynamics reported for high excitation energies in Ref. 26, which are also the dynamics accessed at high temperatures.

Solvent effects on the LiNC ⇌ LiCN isomerization reaction have previously been addressed through the introduction of an argon bath.^16,17,25 The interaction with the bath can be reduced to the Langevin equation through a mapping to the characteristic friction and random noise.^16,19,25 With the approximation $(d/dt)(pϑ/(m*R))≈ṗϑ/(m*R)$⁠, the equations of motion follow from Hamiltonian (2) as

where the stochastic forces ξ_ϑ and ξ_R satisfy the respective fluctuation–dissipation theorems²⁷ 

and

for uniformly distributed noise. The canonical momentum p_ϑ is an angular momentum. It should, therefore, not be surprising that the last term in Eq. (4b) includes a product with the radial coordinate R as it leads to the correct units. The random forces ξ_i are generated at the beginning of the calculation and use a fixed R (= 4.219 6 a.u.) at the barrier with the reduced mass m* (= 2406 a.u.). This is a nontrivial approximation because the reduced mass varies as much as 25% across the positions in R, but it is consistent with prior work, and the error is least when the trajectories are near the barrier. A numerical consequence of this approximation is that the random forces do not vary with R, thus allowing for a significant simplification in coding the equations of motion and in implementing the theory. Thus, the price of this approximation is that the results may be affected in so far as the effective temperature is renormalized.

The results have been calculated with a fourth order Langevin Runge–Kutta method, keeping the random force constant during each time step.

The PGH theory for activated processes driven by Markovian forces has been seen to be very effective in the low- and intermediate-temperature regimes across the Kramers turnover of the rates with respect to friction.²³ Among several examples,^23,28–32 it was shown to be effective for characterizing the dynamics across the potential model of Straub, Borkovec, and Berne (SBB).^33,34 The Langevin implementation is memoryless, which is formally described by the friction kernel,

The SBB approximation uses parabolic functions attached to each other to create a single potential well followed by a saddle as an inverse parabolic function and implemented in the generalized Langevin equation (GLE) with a memory friction. It was used previously by some of us to model the minimum energy path of the LiNC ⇌ LiCN isomerization reaction.²⁵ In the SBB model, memory is introduced through a single exponentially decaying term in the friction kernel,

The propagation of particles inside this potential is strongly dependent on the memory time scale τ = αγ used in the friction kernel and the form of the friction kernel γ(t, t′) itself. For those relaxation processes that occur at times much longer than τ, the response looks ohmic as in the Langevin case. However, the SBB model now allows for dynamical responses from the solvent that can compete with the dynamics in the system. This leads to an effective friction, which arises from the mean of the modulate frictions from the previous times. We found earlier¹⁹ that, in practice, this led to an increase of observed rates by about a factor of 5 when using the Langevin equation (LE) rather than the SBB model.

The first-passage time is the time a particle needs to reach a certain region for the first time, given some initial state. In the case of a reaction, the first-passage time is defined as the time the particle is propagated from a point in the reactant region to a point on some characteristic surface at or beyond a dividing surface (DS).^1–3,5,6,35 In reactive systems characterized by one-dimensional barriers, the DS reduces to a point. It is naively taken to be the saddle point, but other choices are available, just like for variational TST, for example.^36–43 The first-passage times for a series of trajectories from different initial points in the reactive regime experiencing varying thermal forces vary stochastically. Averaged together, they lead to the mean first-passage time t_MFPT, whose inverse is the rate of escape,^6,32,44,45

In the limit of a harmonic barrier, k_MFPT has been seen to be precisely equal to the TST rate, and both are equal to the correct Kramers rate k.³ 

Generally, the rate problem is treated exclusively in the activated regime. Therein, the typical energies of the system are characteristic of an average temperature that is well below the energetic barrier. The smooth turnover in the Kramers rates with friction was resolved by Mel’nikov and Meshkov⁴⁶ and Pollak, Grabert, and Hänggi (PGH).²³ They found a mathematical expression for the rate connecting the low-friction regime—where the rate increases with γ—and the high-friction regime—where the rate decreases with 1/γ. PGH imposes a rate-determining region (or DS in phase space), which requires temperatures to be low enough that the reacting system is somehow thermalized. Since our initial work^17,19,25 demonstrating the applicability of the PGH theory at surprisingly higher temperatures, Pollak and Ianconescu extended it for temperatures near threshold.⁴⁷ 

In the present problem, however, we must also consider much higher temperatures in which the reactive system usually accesses energies much higher than the barrier along the reaction coordinate. The process is consequently effectively barrierless, and the rate problem reduces to the determination of a steady-state flux for a given thermal molecular beam. That is, the typical energies of those states accessing and crossing the barrier at high temperatures—in the sense that they are much larger than the barrier—correspond to states that cross the barrier freely (or ballistically) at their typical velocity. As a consequence, the rates reduce to a simple power law at high temperature,

This behavior is well known since the early work on molecular beams.^48,49 Alternatively, we can recover Eq. (9) from the mean-squared displacement observed in small and long times from the Langevin equation,⁵⁰ where the mean-squared displacement is given by

Thus, the velocity dependence of Eq. (9) is again recovered when

is satisfied,⁵⁰ and the reactants have a short MFPT (and associated small t) due to the fast barrierless crossing.

The PES of LiNC ⇌ LiCN (cf. Fig. 2) is replotted in Fig. 3 in terms of the Cartesian coordinates, x and z, of Li relative to the C_∞v axis of CN. The minimum energy path of the potential nearly follows a semi-circle with a radius of R ≈ 4.5 a.u. The minimum energy of the LiCN molecule (E_min,LiCN = −0.234 21 a.u.) is not as deep as that of LiNC (E_min,LiNC = −0.244 61 a.u.). Consequently, a trajectory starting at the LiCN state has to overcome a smaller barrier height of E^‡ = 0.005 35 a.u. (corresponds to E^‡/k_B = 1690 K) than one starting from the LiNC state. At room temperature (T = 300 K), for example, the barrier of this backward reaction is, thus, low enough that such trajectories are activated to above threshold energy frequently enough that they can be observed numerically. This is also true even for the forward reaction despite its higher barrier.

To ensure that selected trajectories are properly identified as reactive, they must first cross the saddle and then reach the product side without turning around. This condition is satisfied using an absorbing boundary^32,51 defined by an angle of φ = 0.6π and shown as the white dotted line in the Cartesian PES plot of Fig. 3. Indeed, if the trajectory reaches this line, then it generally has enough momentum in the direction of the product side to make it very unlikely for it to turn around and climb back over the saddle to the reactant side. Use of absorbing boundary conditions along space and/or energy to address low to mid friction also has significant precedent in the literature.^23,46,52 In the present case focusing on high temperatures, we did not add an energy condition to the absorbing boundary because the spatial constraint is sufficient to remove the high-energy trajectories once they have reached the product basin (cf. Fig. 3) before their long-time return to reactants. Its use enables us to better compare the MFPT calculations to the AAMD results obtained earlier in Ref. 16 where direct rates were determined based on the flux through a DS along the reaction coordinate.²⁵ 

The representative trajectory shown in Fig. 3 is only one of the 1500 propagated trajectories from the ensemble used to calculate the MFPT rate k_MFPT. Trajectories of this thermal ensemble are initialized at a specific temperature T and located near the reactive well. For simplicity in the current implementation, they are all placed at the LiCN minimum at R = 4.7947 a.u. and ϑ = 0, which in Cartesian coordinates corresponds to x = 0 a.u. and z = 4.7947 a.u. Their velocity is sampled from a Maxwell–Boltzmann distribution with temperature T. The assumption here is that the redistribution to a quasi-bound distribution in space is fast compared to the escape times when the system is trapped. Meanwhile, in the barrierless cases, there would be no such trapping and we are merely considering the steady-state flux from the injection of particles at this origin. This case could presumably be accessible through increasingly exquisite spectroscopic methods, such as transition state (TS) spectroscopy.^53–56 

At higher bath temperatures, the distributions in energy—total, potential, and kinetic—of the initial reactant ensemble necessarily shift to higher values. Consequently, a larger percentage of the states in the reactant ensemble will possess the required energy to surmount the barrier at any given time. More of these states escape over the barrier and are removed from the simulation. The remaining and decreasing population of reactive states is re-thermalized according to the bath temperature and thereby continues to contribute to the steady-state rate. As the temperature further increases, this activation is faster, leading to increasingly larger kinetic energies in the trajectories contributing to the steady-state rate. This, in turn, causes smaller MFPTs and higher rates. Specifically, ensembles thermalized at some temperature T and propagated at a certain friction value γ lead to rates that depend not on the fast inertial trajectories but rather those that endure several traversals of the reaction region—at energies below activation threshold—before escaping.

We now report the MFPT rates k_MFPT across the friction domain at several intermediate temperatures in Fig. 4. As before,^16,19,25 we observe a clear Kramers turnover, and the shape of k_MFPT is in very good agreement with the corresponding PGH rate formula. The rate maxima are always around γ ≈ 4 × 10⁻⁴ a.u.

To obtain the rates in Fig. 4, different ensembles at T = 300 K, T = 450 K, and T = 600 K are thermalized on each friction value γ. Each mean rate k_MFPT is calculated using 1500 trajectories propagated by the Langevin equation. At low friction (1 × 10⁻⁵ a.u. ≪ γ ≪ 2 × 10⁻⁴ a.u.), we find the expected linear increasing behavior of the rates. Similarly, at high friction, the rates decrease strongly with 1/γ. These two limits have been known since the work of Kramers. In combination, they give rise to the eponymous Kramers turnover. In absolute values, however, the MFPT rates overestimate the expected rates by a factor of 5. This is unfortunate, but not unexpected, because the MFPT rates—just like those from flux over population³—are known to be approximate due to several factors. In the present implementation, these MFPTs are uniformly shorter, leading to faster rates than measured through the other approaches, presumably because there is an acceleration from the initial distribution that is not being dissipated by the bath.

The rates calculated using trajectories propagated via the LE are also overestimated by about a factor of 5 compared to the SBB model and AAMD calculations. Specifically, the Kramers turnover rate maxima k_max in the LE and generalized Langevin equation (GLE) from Ref. 19 are in agreement only if the rates obtained for the reactive flux calculations for the LE were to be divided by a factor of 5. Meanwhile, the GLE rates obtained using PGH theory for the near ohmic limit (at very small α) do not suffer from this overestimation. This suggests that the LE rates calculated using reactive flux suffer from numerical error not present in the GLE dynamics. We conjecture, though have not proven, that this discrepancy arises from a slower energetic redistribution among the quasi-bound states, allowing for a faster escape of the initial higher energy configurations.

To confirm the effect of friction in the Kramers turnover of k_MFPT seen in the PGH theory in Fig. 4, the bath parameter α in the PGH friction kernel is set to a small value α = 1.5625 a.u. The maximal rates for fixed α in this case are still found at friction values γ₀ in a similar range, and hence, this case (α = 1.5625 a.u.), indeed, leads to a short memory time scale of τ ≈ 6.25 × 10⁻⁴ a.u. This places it in the memoryless dynamics regime characteristic of the LE, as noted above. Meanwhile, the PGH rates obtained at the longer memory time τ ≈ 5.25 × 10² a.u. —corresponding to a larger α—as was observed in the AAMD simulations, still exhibit similar maximal rates in Fig. 5. At intermediate temperatures, for example, the memory time scale is still much lower than the MFPT (t_MFPT ≳ 1 × 10⁴ a.u.), and thus, the system remains in a nearly memoryless regime. As the PGH theory—which is also known to be in good agreement with TST³² at the appropriate limits—describes the low- and intermediate-temperature regime well, it is not surprising that the MFPT rates are effective as long as we ignore a temperature-independent uniform factor.

The MFPT rates can also be used to describe the high-temperature regime, where the typical energies are well above the barrier height E^‡/k_B = 1690 K. Therein, the Kramers turnover maxima in the MFPT rates are compared to the corresponding reactive-flux rates¹⁶ in Fig. 5. Specifically, we compare against LE-, GLE-, and AAMD-based calculations, where we consider AAMD to yield the most accurate results. Notably, the latter includes cavity reorganization effects, which may not be fully described within the LE-based reactive-flux and MFPT approaches.

Figure 5 also includes TST benchmark calculations based on the Eyring–Polanyi (cf. Appendix B 2) and the Polanyi–Wigner (cf. Appendix B 1) rate expressions. While the Eyring–Polanyi rate describes the reaction rate reasonably well for k_BT ≪ E^‡, as is well known,^{32,39,57–63} it quickly deviates from the other results as the temperature is increased. Indeed, the usual rule of thumb is that there is a breakdown in rate theories—because they begin to violate the separation between reactant and transition state regions—above this threshold. Thus, the result in the all-atom LiCN dynamics in an argon solvent seen in Refs. 16 and 25 was surprising because it suggested the existence of an observable rate at such high temperatures, and their turnover with an apparent friction obeyed the Kramers turnover seen through the lens of PGH theory. Here, we see that part of the origin of the breakdown in the Eyring–Polanyi theory is that it suggests an infinite rate, which is not possible when the particles move at velocities distributed around a Boltzmann distribution for a given k_BT. Instead, at high enough temperatures, the reactants move at apparent rates that are proportional to their average velocity because at such energies, the barrier is not visible. The results in Fig. 5 demonstrate that the rates, as also captured by the Polanyi–Wigner expression, move smoothly from the well-known low-temperature regime through intermediate temperatures in which the barrier plays a role, albeit a decreasing one with increasing temperature.

Thus, at increasing temperatures, the reaction is no longer limited by the activation process that is inherent in the Eyring–Polanyi formula and rises much faster, leading to a very large upper bound of the rate. On the other hand, the Polanyi–Wigner formalism appears to capture the slowing down in the rates due to the decreasing quasi-bound population in the reactive basin with increasing temperature. Specifically, the prefactor in the Polanyi–Wigner rate is no longer temperature dependent but rather determined by the vibrational frequency ω₀/2π (⁠$≈$3.2 × 10¹² s⁻¹) of the reactant well and, therefore, leads to a flattening in the rate with increasing high temperatures. Nevertheless, the interaction with the solvent manifests itself by way of the decreased rates through the intermediate temperature regime.

The comparison between the maxima of the PGH rates and the AAMD rates for higher temperatures in Ref. 16 is also shown in Fig. 5. The agreement shows that the PGH theory can also be applied to high temperatures. In this work, we further found that the PGH rates are able to describe the intermediate- and high-temperature regime, both with the memory time scale appropriate to the system (cyan pluses) and through the approximate Langevin approach (purple crosses). This also justifies the use of the k_MFPT obtained from the LE friction kernel at the high temperatures.

The AAMD, GLE reactive-flux, and PGH results reported in Fig. 5 are in good agreement at high temperatures. As already observed in Ref. 19, the reactive-flux results based on the LE grossly overestimate the rates by a factor of 4.5–4.8 and 5.1–9.3 compared to those seen in the GLE and AAMD reactive-flux calculations, respectively. The MFPT rates for the LE are lower than those from the reactive-flux calculations on the LE but not enough to match the other methods. Specifically, k_MFPT overestimates by a factor between 2.9 at T = 2500 K and 2.2 at T = 5500 K compared to the GLE reactive flux. The overestimation factor is larger when compared to AAMD calculations, monotonically rising from 3.3 at T = 2500 K to 4.6 at T = 5500 K. There is, however, no such simple trend when considering the whole temperature range via comparison with PGH rates. Instead, the overestimation factor is particularly high at very low temperatures (9.7 at T = 300 K), quickly drops until it reaches a local minimum at T = 2500 K, and finally rises again slowly with temperature.

Compared to the LE reactive-flux approach, the MFPT appears to better capture the escape of the quasi-trapped initial states at high temperatures. These low-probability states contribute to the steady-state flux obtained from the AAMD, GLE reactive flux, and PGH, leading to an even lower rate. The disagreement between the latter and those obtained using the LE is surprising because the PGH results for a GLE with a very small α are formally the exact answer for an LE and should be in direct agreement with the reactive flux obtained on the LE. Thus, both the MFPT and reactive-flux calculations of the LE may be suffering from numerical errors due to the fact that the needed contributing trajectories are rare. The better agreement obtained with the MFPT may reflect the fact that it does not suffer as much from the constraints related to the selection of the initial states from the quasi-bound population.

Despite some of the loss in accuracy arising from the strong localization of the initial distribution of states, the advantage of the MFPT rates compared to AAMD rates is that they can be calculated for the whole temperature regime in less time. Furthermore, a square root function shape is observable for k_MFPT, as described in Sec. II D. This shape arises because the rate is known to depend on the square root of temperature in this regime.

The difference between k_MFPT and reactive-flux-based approaches in Fig. 5 is a consequence of the assumptions inherent in each method. For the reactive-flux rate, an extended ensemble is initialized at the minimum of LiCN and the rates are obtained by the flux-over-population method across numerically integrated trajectories. In the GLE, the rates are obtained through a direct (and analytic) analysis of the reactive ensemble across the barrier. In principle, long-time returns to the barrier may be absent from the latter, leading to an overestimate in the transition state theory calculation. Such long-time sojourns are long only in the sense that they are longer than the effective reactive timescales but may be short in an absolute sense because the temperatures (and the kinetic energies) are large. A related consequence of the high-temperature regime is that the kinetic energy tends to be dominated by equipartition, leading the effective average rate across the barrier to simply be described by the barrierless average velocity. The ensuing square root behavior is seen correctly in the MFPT rates in Fig. 5 at high temperatures, whereas the PGH rates are unfortunately linear.

In this work, we have calculated the rates of the LiCN → LiNC backward reaction using mean first-passage times. We benchmarked them relative to the results of earlier AAMD calculations and the PGH rates across several temperatures and friction regimes. The MFPT rates are effective in describing the Kramers turnover in the intermediate temperature regime as shown by comparisons to the PGH rates calculated with a short memory time scale. At very high temperatures, when the reactants motion effectively becomes ballistic, the MFPT approach was seen to correctly capture the square-root dependence of the rates on temperature.

At low temperatures, there are several examples of near-perfect unadjusted agreement between methods such as the reactive flux and GLE-based theories such as Grote–Hynes theory.^64–66 See, for example, reports of the rates in cyclohexane interconversion,⁶⁷ NaCl dissociation,^68,69 S_N2 reactions,^70–72 and even some enzymes.^73–75 However, for the LiCN isomerization reaction at intermediate and high temperatures, we found, as before,¹⁹ that methods based on Langevin dynamics overestimate the rate roughly by a factor of 5 compared to AAMD or GLE simulations. The program used for our paper was written independently of the one used in Ref. 19. It, therefore, seems unlikely that this factor is merely caused by a bug in some code. Instead, this suggests a physical origin which we conjecture here, as done earlier in Ref. 19, that it originates in otherwise unaccounted mean field dissipation in the determination of the stochastic trajectories. Thus, the MFPT approach pursued in this paper was found to provide no worse agreement than other methods while being relatively simple to implement in systems with arbitrary dimensionality. This overestimation, nevertheless, poses a challenge for future work.

Furthermore, comparisons of the computed rate constants to an experiment would clarify which theoretical approach is most suitable to describe the LiCN isomerization reaction. The experimental measurement of these rates, especially in the high-temperature regime, also poses a challenge for future work. Specifically, for LiCN, it may be difficult to construct an initial distribution of states localized at the reactant well that is weakly solvated by a thermalized argon bath. Perhaps, spectroscopic scalpels such as those offered by TS spectroscopy^53–56 could be employed, and we look forward to seeing such advances.

In addition, the MFPT approach provides a deeper view of the structure of the LiCN isomerization reaction. The footprint (or projection) of the ensemble of trajectories that contribute to the rate onto the domain of reactant coordinates appears to spread across the intrinsic reaction coordinate (IRC).^76,77 Although not addressed quantitatively here, we conjecture that this spread can be used to characterize the manifold (or tube) that envelopes the possible reaction pathway⁴⁰ but now constructed organically without the assumptions required of the IRC.

This work also opens at least two new possible directions to address rates in chemical reactions and other activated systems: The first lies in the need for a more direct comparison between NHIM rates²⁰ with MFPT rates. The second lies in the use of these approaches to obtain the rates for a time-dependent LiCN system driven by an electrical field so as to fully demonstrate the computational advantages of this approach in chemical reactions generally. Finally, reproducing the calculations of the MFPT on the Langevin equation with those for the general Langevin equation (with the observed memory time of the solvent) could also confirm the generality of our finding.

Useful discussions with Thomas Bartsch are acknowledged. The German portion of this collaborative work was partially supported by the Deutsche Forschungsgemeinschaft (DFG) through Grant No. MA1639/14-1. The Spanish portion was supported by the Spanish Ministry of Science, Innovation and Universities (MICIU) under Contract No. PGC2018-093854-B-I00. The U.S. portion was partially supported by the National Science Foundation (NSF) through Grant No. CHE 1700749. This collaboration also benefited from the European Union’s Horizon 2020 Research and Innovation Program under Marie Skłodowska-Curie Grant Agreement No. 734557.

The authors have no conflicts to disclose.

M.M.S. and J.R. contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The LiNC ⇌ LiCN PES according to Ref. 21 consists of two parts: a damped long-range energy plus a short-range energy.

The long-range part is composed of the electrostatic energy

and the induction energy

Here, P_L is the Legendre polynomial of order L, ⟨Q_L,0⟩ denotes the expectation value of the order-L CN⁻ multipole moment, and $Cl1,l2,L$ are the induction coefficients. Numerical values for ⟨Q_L,0⟩ and $Cl1,l2,L$ are given in Table I. The damping is represented by

with fit parameters a = 1.5156 a.u. and R₀ = 1.9008 a.u.

The short-range term can be written as

where D_L(R) has been fitted to the analytical form

Numerical values for the fit parameters A_L, B_L, and C_L are given in Table II.

Combining long and short-range energies, the final PES reads

We noted a discrepancy while comparing the minimum energy path calculated using the parameters published in Ref. 21 with Fig. 2 from Ref. 21. As shown in Fig. 6, the two curves differ visibly. Two parameters in the original source code⁷⁹ differ significantly from the originally published values, possibly due to errors introduced during the paper’s production process. The correct values—shown in bold in Tables I and II—yield a much better agreement with Fig. 2 from Ref. 21.

An implementation of the potential in the Python programming language using the same parameters as the original source code⁷⁹ can be found on GitHub.⁸⁰ 

Most approximative rate formulas follow the Arrhenius form^32,81

where ν is a possibly temperature-dependent prefactor, E^‡ is the reaction’s barrier height or activation energy, k_B is the Boltzmann constant, and T is the temperature. In the following, we present two important TST variants of this equation.

One of the earliest results of TST is the unimolecular rate equation derived, among others, by Polanyi and Wigner in 1928.⁸² It follows the Arrhenius rate law (B1) with the pre-exponential factor ν(T) given by the vibrational frequency ω₀/2π of the reactant well,³² 

where E^‡ is the internal energy difference between the barrier and the reactant of the isolated system. This equation has been derived in various contexts. It can, e.g., be recovered from the underdamped regime ω^‡k_BT/E^‡ ≪ γ ≪ ω^‡ of Kramers’s medium-to-high-viscosity rate,^24,46

where γ is the friction and ω^‡ is the inverse barrier frequency. Equation (B2) can, therefore, be seen as a TST rate, which is an upper bound for the rate at the turnover in Kramers’s original theory for a solvated reaction.³² 

The usual (or modern) form for the classical TST rate equation is given by the Eyring–Polanyi equation^36,83,84

with κ = 1, where κ is the transmission coefficient, h is Planck’s constant, and ΔG^‡ is the Gibbs energy of activation in the context of a solvent.

The Gibbs energy of activation can be approximately determined from the enthalpy of activation ΔH^‡ via

where ΔS^‡ is the entropy of activation. In turn, the enthalpy of activation for an unimolecular gas-phase reaction can be written as

The energy of activation E^‡ and the entropy of activation ΔS^‡ can finally be determined from the minimum energy path and the potential of mean force of the reaction^16,25 by equating the latter with the Gibbs energy.

The transmission coefficient κ in the Eyring–Polanyi equation describes the fraction of states that cross the DS between reactants and products at most once, i.e., those that do not recross. This quantity cannot be determined from straightforward statistical mechanics and is, therefore, of great interest⁸³ in the general case. In TST, it is assumed to be ∼1.³² This approximation is valid if the temperature is not too high or if the friction is sufficiently strong.