Competing quantum effects in heavy-atom tunnelling through conical intersections

Thermally activated chemical reactions are typically understood in terms of overcoming potential-energy barriers. However, standard rate theories break down in the presence of a conical intersection (CI) because these processes are inherently nonadiabatic, invalidating the Born–Oppenheimer approximation. Moreover, CIs give rise to intricate nuclear quantum effects such as tunnelling and the geometric phase, which are neglected by standard trajectory-based simulations and remain largely unexplored in complex molecular systems. We present new semiclassical transition-state theories based on an extension of golden-rule instanton theory to describe nonadiabatic tunnelling through CIs and thus provide an intuitive picture for the reaction mechanism. We apply the method in conjunction with first-principles electronic-structure calculations to the electron transfer in the bis(methylene)-adamantyl cation. Our study reveals a strong competition between heavy-atom tunnelling and geometric-phase effects.

We add an additional mode, xb , that couples to the tuning mode with coupling strength, c b , and approximately accounts for the frictional effect of the remaining molecular modes.The potentials as functions of these three modes are defined as where the positive (negative) sign corresponds to the reactant (product) surface and the parameters are defined in Table SI.1.The parameters are mass-weighted such that m is arbitrary and can thus be chosen as m = 1.The diabatic coupling is ∆(x) = α • x = αx c , where α = α .

SI.2 Rate theory
The diabatic Hamiltonian for a general nonadiabatic reaction is given by The nuclear Hamiltonian of electronic state n ∈ {0, 1} is of the form Ĥn = f j=1 p2 j 2m + V n ( x), (SI.3)where x = (x t , xc , x3 , . . ., xf ) is the coordinate vector comprising the tuning mode (x t ), the coupling mode (x c ) as well as f − 2 further nuclear degrees of freedom.Without loss of generality, the nuclear degrees of freedom have been mass-weighted such that each has the same mass, m.The nuclei move on the potentials V n ( x) with conjugate momenta pj .

SI.2.1 Quantum-mechanical and classical rates
The full quantum-mechanical thermal reaction rate constant is given by [3][4][5] k Z 0 = For the simple harmonic model [Eq.(SI.1)], the exact rate can be computed by expanding the Hamiltonian and flux operators in a basis of vibronic states.The linearly coupled systembath model in Eq. (SI.1) can be transformed to a normal-mode basis which diagonalizes the potential energy.In the basis of normal modes, {q i }, with corresponding frequencies {Ω i } and shifts {ξ i }, the potentials are given by where again the positive (negative) sign corresponds to the reactant (product) surface.Note that because the coupling mode of the system-bath model is not coupled to the other modes, it will not be altered by the normal-mode transform and thus ω c appears in the set {Ω i }.For both correlation functions, we choose τ = β /2.To obtain a well-defined rate constant that can be compared to the instanton results, we integrate the correlation function numerically over the interval [−t p , t p ], which encompasses the central peak only.Here we choose a plateau time of t p = 0.05 ps.We thus ignore resonances at later times, which we assume will be damped out by friction from the other modes in the fulldimensional system.For the calculation of the full quantum-mechanical correlation function, we employed a basis of 9 × 15 × 30 vibrational wavefunctions for the tuning, coupling and bath modes on each electronic state.
Since the three normal modes are uncoupled, the vibrational eigenstates (ψ µ 0 , ψ ν 1 ) of the nuclear Hamiltonians, Ĥn , with corresponding energies (E µ 0 , E ν 1 ) can be constructed as direct products of the eigenstates of the individual modes.The vibronic states are formed as a direct product of the vibrational states with electronic states |n .
We then evaluate the matrix exponentials in Eq. (SI.4) in a truncated vibronic basis, enabling the calculation of the flux correlation function.The result needs to be converged with respect to the size of the vibronic basis.The final time integration over the correlation function is carried out numerically.The integration limits are reduced to [−t p , t p ] so as to encompass the first plateau only, as illustrated in Fig. SI.1.This permits the definition of a rate in a finite-dimensional bound system.

SI.2.1.1 Fermi's golden rule
We now consider the golden-rule approximation.For small coupling strength, ∆( x), Eq. (SI.4) can be expanded in a perturbation series and approximated by the first-order term, which results in the following expression for the golden-rule rate constant 6,7 with the golden-rule flux correlation function, c GR (τ + it).The trace, tr{•}, is over nuclear degrees of freedom only.To arrive at the famous Fermi's golden-rule rate 8,9 generalized for an initial thermal distribution, 10,11 one evaluates the trace in the basis of reactant and product vibrational eigenstates and takes the time integral analytically In order to define a rate in the bound three-mode model defined above, one can employ the function , (SI.9) which effectively smears the harmonic-oscillator energy levels.

SI.2.1.2 Path-integral formulation of Fermi's golden rule
Instead of using the standard Fermi's golden rule, we wish to find an expression for the rate that does not require knowledge of the vibrational wavefunctions.Inserting the linear form of the diabatic coupling at a CI into the correlation function in Eq. (SI.7) leads to (SI.10) The quantum golden-rule rate for our harmonic model system with linear diabatic coupling can equivalently be evaluated using path-integral theory.The expression for the pathintegral propagators for a system with the potentials in Eq. (SI.6) are known 12,13 and the position integrals with monomial prefactors in Eq. (SI.10) can be solved analytically, leading to the rate expression where For the charge transfer in the BMA cation, the thermodynamic driving force, ε, is zero and the value of τ that makes φ(τ ) stationary is τ = β /2 by symmetry.However, we have presented the equations such that they are valid also for systems where ε is non-zero.
We illustrate the GR flux correlation function in Fig. SI.1, which is seen to be in excellent agreement with the non-perturbative correlation function.In order to obtain a well-defined rate constant, we again integrate over the interval [−t p , t p ].

SI.2.1.3 Classical golden rule
The classical golden-rule rate can be computed by evaluating the trace in Eq. (SI.10) by a phase-space integral where p is the vector of momenta conjugate to the coordinates x.The Gaussian integrals over the momenta are easily solved analytically and the integral over time results in a δ-function.
In general, the position integral cannot be solved analytically and is thus approximated by the method of steepest descent.Thereby, one locates the stationary point of the effective potential V = (1 − λ)V 0 + λV 1 in the combined space of positions and the Lagrange multiplier λ = τ /β . 14The stationary point of V corresponds to the minimum-energy crossing point (MECP) of the two potentials or, in our case, the minimum-energy conical intersection (MECI).It can be viewed as the golden-rule analogue to a transition state in adiabatic reactions.Carrying out the position integrals leads to the following expression for the classical golden-rule rate 13,15 where V ‡ is the potential at the MECI relative to the potential of the reactant minimum, For globally harmonic model systems such as in Eq. (SI.1), the classical rate [Eq.(SI.14)] simplifies to where the Marcus-theory Franck-Condon factors f FC,MT are given in terms of the reorgani- , and the thermodynamic driving force.In this way the standard Marcus-theory rate (for constant non-zero ∆) given by is generalized to describe the effect of the CI using the prefactor χ MT . 16,17

SI.2.2 Instanton rates SI.2.2.1 Semiclassical instanton theory
The SCI rate for a reaction proceeding via a CI is derived as a semiclassical approximation to Eq. (SI.10).We evaluate the trace in the position representation, insert a complete set of position states in the middle and replace the resulting quantum propagators by semiclassical van-Vleck propagators, which leads to where the prefactors, C n , are assumed to vary slowly compared to the exponential and can thus be pulled outside of the integral.In order to carry out these integrals by steepest descent, we locate the stationary point of S (called the instanton) in the combined space of positions and imaginary time and keep t = 0 since this is the dominant point of the correlation function.The action can then be expanded in a Taylor series to second order around the stationary point, which gives where x = (x , x ) is a vector combining all 2f position coordinates and ∂ 2 S ∂x∂x is the secondderivative matrix of the action with respect to these coordinates evaluated at the stationary point with positions x = ( x , x ) and imaginary time τ .What remains is a standard multidimensional Gaussian integral with monomial prefactors, which evaluates to where C is the determinant of the second-derivative matrix of the action, C = ∂ 2 S ∂x∂x .The two-point correlation χ arises from the monomial prefactors inside the integral and is given by χ = (α ) T C −1 α , where α = (α, 0) and α = (0, α) are vectors in the space of end points (x , x ) and 0 is an f -dimensional vector of zeros.The remaining time integral over the correlation function can also be evaluated by steepest descent, which is equivalent to the SCI approach without a CI and again assumes the prefactors, this time including χ, to vary slowly in time compared to e −S/ . 13The full SCI rate constant is then given by where we used the Cauchy-Riemann equations

SI.2.2.2 Steepest-descent instanton theory
The calculation of the SDI rate is conceptually similar to the calculation of the standard golden-rule instanton rate without CI with the exception that the diabatic coupling is included in the action [Eq.( 4 The total rate is thus given by a sum over the steepest-descent rates associated with each of these stationary paths e −S eff ( x a , x a ,τa)/ , (SI.20)   where all quantities are evaluated for the respective stationary path of S eff .Here, the C eff prefactor is defined in analogy to C in Eq. (SI.18) but with respect to the effective action in- To construct the plot, we transformed the system to a spin-boson model, as described in Eq. (SI.6), and used the corresponding analytic expressions for the action. 13he two remaining coordinates, as well as the value of τ were fixed to their values at the stationary point of the (effective) action.Note that in this symmetric spin-boson model, the stationary points in the SCI and SDI formalisms only differ in the coupling coordinate, while the values of the two other modes and τ are the same.stead of S. We thereby separated the range of the position integrals into different regions, Γ a around the stationary points (similarly to Ref. 18).The prefactors, χ eff a = sgn[∆( x )∆( x )], depend on the relative sign of the diabatic-coupling terms at the hopping points of the stationary paths.In particular, direct instantons for which x = x have χ eff a = 1, whereas those which wind around the CI have χ eff a = −1.The derivatives of the coupling also enter into the second-derivative matrix C eff and the total τ -derivative since these terms derive from S eff .
In Table 1 of the main text, we presented the SDI rate constant for the three-mode model defined above.Note that the 20% deviation from the quantum-mechanical result arises from the neglect of anharmonicity around the stationary points of S eff , which within SDI is not quadratic even for the harmonic LVC model, because of the contribution from the diabatic coupling.

SI.2.3 Resummed semiclassical instanton theory
We presented the SCI result for the three-mode model in Table 1 of the main text, which deviates from the FGR rate constant by a factor of around 2. Considering all error sources in realistic applications, including the electronic-structure calculations, this still constitutes decent agreement.However, we can understand the reason for the deviation by comparing the flux correlation functions underlying the different methods, which are illustrated in Based on our understanding of the error in SCI, we can improve the method by at least partially correcting for a rapidly oscillating prefactor.First, we go beyond the constant approximation for χ by expanding it as a second-order Taylor series.Next, we can use our knowledge that the analytical expression for χ in a harmonic system [Eqs.(SI.12a)] is a cosh function.Given χ(τ ) and its second τ -derivative we can therefore propose a resummation for the Taylor series as where χ = ∂ 2 χ ∂τ 2 , and we used the relation cosh(it) = cos(t).The corresponding approximation of the correlation function, depicted in In order to carry out this approach for an anharmonic system, we need to compute the second time-derivative of χ at τ = τ .This can be achieved by taking time-derivatives of the correlation function, as can be seen from Eq. (SI.18).Taking the second time-derivative of the golden-rule flux correlation function [Eq.(SI.10)] leads to where the second term on the right-hand side arises from time-derivatives of the timeevolution operators and is thus also found in Wolynes theory, 6,20,21 even in cases where the position-dependence of the diabatic coupling is ignored.It follows that the first term on the right-hand side of Eq. (SI.22), which contains time-derivatives of the diabatic coupling, contains the information about χ.In particular, the term of interest is given by As before, we expand this expression in a position basis, replace the quantum propagators by their semiclassical analogues and evaluate the term around the instanton (the stationary point of the total action).At the instanton, the component of the momentum pointing along the coupling mode (orthogonal to the reaction coordinate given by the tuning mode) is zero.
We therefore need to use the following expansion where all derivatives are evaluated at the hopping points.The analogous procedure is followed for the terms associated with the propagation in the other direction.Second derivatives of S n can be evaluated using formulas presented in the Appendix of Ref. 22.These require knowledge of nothing more than the potentials, gradients and Hessians along the instanton path, which are anyway required for the original SCI approach.
We are now in the position to evaluate χ at the instanton, which is given by χ and the masses are accounted for by mass-weighting all terms.With the general expressions for χ(τ ) and χ(τ ) at hand, the resummed SCI formula for the rate constant is given by where χ r = χ e 2 χ χ S−1 , and all quantities are evaluated at the instanton.We present the rate constant for the three-mode model obtained with this method in Table SI.2.It is seen to be in excellent agreement with the quantum-mechanical results, as expected for the harmonic system.For more complicated, anharmonic systems where the diagonal and off-diagonal blocks of the C-matrix cannot be simultaneously diagonalized, the resummation would have to be adapted.We leave the development of such a more general resummation open for future work since the correction of the SCI rate in the model system is relatively minor.

SI.2.4 Extended results for the model system
In Table SI.2, we compare the classical and NA-TST rates to the other rates given in Table 1 of the main text.Furthermore, we illustrate the rate constants computed with several methods over a range of different temperatures in Fig. SI.4.As also seen in our ab initio results, the NA-TST rate for the model system can be larger than the instanton and quantummechanical rates, despite the neglect of nuclear tunnelling in this method.Furthermore,

SI.3 SDI and the geometric phase
It is also possible to define a classical rate in the spirit of the SDI approach.Therefore, we evaluate the correlation function [Eq.(SI.7)] by a phase-space integral but this time include the coupling in the exponent.We then arrive at a configuration-space integral that we can take by steepest descent about the stationary points of the effective potential /β.This leads to the rate expression where the sum is over stationary points of the effective potential, the MECI partition function, Z eff, ‡ , is evaluated based on the second-derivative matrix of V eff with respect to positions only with the reaction coordinate projected out.Finally, χ eff cl = 1 is the classical limit of χ eff .For the model system from Sec. SI.1, the method leads to a rate of 1.010 ns −1 , which is in excellent agreement with the classical rate in Table SI.2 obtained with Eq. (SI.14).Note that in the effective SDI action [Eq.( 4)], ∆(x ) and ∆(x ) can vary independently giving rise to two pairs of stationary paths, one pair with a positive sign and the other pair with a negative sign.In contrast, the ∆(x) 2 term in the effective potential of Eqs.(SI.27)only allows for one pair of stationary paths with a positive sign.The classical rate therefore does not contain negative contributions and cannot account for the GPE.The lack of GPEs in the classical and NA-TST rates is the reason why these rates can be higher than the instanton and quantum-mechanical results, as shown in  Figure SI.5:Contour plot of the adiabatic ground-state potential projected into the x t x cplane (with a relaxed bath-mode coordinate) along with the steepest-descent instantons and the stationary points of V eff (black crosses) defined beneath Eq. (SI.27) at 300 K.The direct blue and red paths contribute to the rate with a positive sign, whereas the winding violet paths (clockwise and counterclockwise, lie on top of each other) contribute with a negative sign.Additionally, the Born-Oppenheimer instantons (black lines) and transition states (black dots) are shown.Note that there are two BOI paths curving left and right around the CI, which are not resolved here because the paths traverse the barrier very close to the CI.
It can further be seen that the curving SDI paths deviate significantly from the corresponding classical stationary points.In both the instanton approach and the classical approach, a compromise needs to be found between maximizing the diabatic coupling, which grows away from the CI, and minimizing the potential-energy barrier, which is smallest at the CI.Such effects of corner cutting, where tunnelling changes the reaction mechanism, are not uncommon.Usually, however, the instanton cuts the corner to find a shorter, more direct path from reactant to product at the expense of travelling through regions of higher potential energy.This can be more favourable because tunnelling through high but narrow barriers is very efficient.In this case, the steepest-descent instantons take a longer path exhibiting a higher barrier than the path via the classical stationary points, which is favourable because of the larger diabatic coupling along these paths.The ability of the instanton to describe tunnelling hence allows it to find a better (in fact the optimal) compromise between path length, potential-energy barrier and diabatic coupling than classical rate theory.

SI.4 Ab initio calculations
In this section we present additional details of the ab initio calculations and additional discussions on the ab initio results.

SI.4.1 PBE vs PBEX50 functional
Hybrid functionals, often with an increased fraction of Fock exchange, are commonly used in constrained density functional theory (CDFT) calculations, especially for computing diabatic coupling. 23,24Since the PBE functional is considerably less expensive than the PBEX50 functional (i.e. a hybrid functional with 50% Fock exchange), we employed it for the instanton optimizations.We have carefully tested its performance compared to the PBEX50 functional.
The fully optimized geometries with the PBEX50 functional are almost identical to those of the PBE functional, and the MECP barrier is 0.180 eV (only 0.009 eV higher than the PBE barrier).Therefore we think that it is justifiable to use the PBE functional in the geometry (including instanton) optimizations.

SI.4.2 Coupling mode from CDFT
Since the coupling mode vector is defined as α = ∂∆ ∂x , it seems one could straightforwardly perform a finite displacement calculation at the MECI to obtain α.However, CDFT is only able to return unsigned |∆| values, making the procedure to obtain α more complicated as |∆| is not a differentiable function at the CI.We first perform forward-difference and backwarddifference calculations, and by taking the unsigned average of the two (i.e., the forwarddifference derivative has only positive elements and the backward-difference derivative has only negative elements so we reverse the backward-difference derivative before averaging), we obtain α = ∂|∆| ∂x .This vector has the correct magnitude (i.e. the same as α), but the relative signs for each element are missing.Next we displace the MECI in an arbitrary direction by a displacement larger than the step size used for the numerical derivative.At this point, we perform a central-difference gradient calculation to obtain the signs of the elements of α.The procedure is repeated to ensure consistency by recomputing the gradient at a point displaced from the CI in the direction of α, where the derivative of |∆| is well defined.Finally, we checked that displacing along α does not lift the degeneracy of the two diabatic states and that the slope of |∆| along this mode is indeed α = ||α||.

SI.4.3 Rate constants and H/D kinetic isotope effects
Here we present all the computed rate constants for BMA and deuterated BMA.One can see from Table SI.3 (comparing the top and bottom sections) that the H/D kinetic isotope effect (KIE) predicted by NA-TST is 1.2 at 300 K, while the KIE predicted by SCI is 1.4.
At 150 K, KIE only slightly increased (to 2.3) according to SCI, only slightly larger than the NA-TST KIE (1.4).By comparing the decompositions in Table SI. 4 and Table 2 in the main text, one can see that apart from the ZPE contributions to the H/D isotope effect (which appear in the NA-TST KIEs), the remainder mainly originates from the tunnelling factor.Meanwhile, if one substitutes all 12 C atoms with 13 C, the tunnelling factor decreases from 9.1 to 8.6 at 300 K, which as discussed in the main text, is a remarkably large KIE for heavy-atom substitution at room temperature.Given the seemingly negligible corner-cutting effects, one might expect that a one-dimensional description based on the minimum-energy pathways (MEPs) would be sufficient.However, the contribution from fluctuations (i.e. the ratio between SCI and NA-TST prefactor) at 150 K would be 0.40 if computed on the one-dimensional MEPs, significantly larger than the value (0.023) computed on the full-dimensional system.One can see from Table SI.5 that the contribution from fluctuations decreases rapidly with temperature for the full-dimensional system, while it changes much slower for the 1D MEPs.This indicates that a theory based on the one-dimensional MEP is inadequate to describe the system, and that the fluctuations around the instanton contain important multidimensional effects.A similar behaviour of the fluctuations has also been found for the LVC model (see Table SI .6).This indicates that this effect is not unique to the ab initio system and that it can be reproduced with multidimensional globally harmonic models.

SI.5 Ab initio system vs LVC model
In  highlights that the tunnelling contribution constitutes the decisive difference between the two approaches.This deviation originates from differing instanton actions, indicating that the barrier in the harmonic LVC model is narrower than the one obtained with on-the-fly CDFT.

Figure SI. 1 :
Figure SI.1:Full (exact) quantum-mechanical and golden-rule (GR) flux correlation functions as a function of time.For both correlation functions, we choose τ = β /2.To obtain a well-defined rate constant that can be compared to the instanton results, we integrate the correlation function numerically over the interval [−t p , t p ], which encompasses the central peak only.Here we choose a plateau time of t p = 0.05 ps.We thus ignore resonances at later times, which we assume will be damped out by friction from the other modes in the fulldimensional system.For the calculation of the full quantum-mechanical correlation function, we employed a basis of 9 × 15 × 30 vibrational wavefunctions for the tuning, coupling and bath modes on each electronic state.

G 0 and G 1
are the gradients on the two potentials and H = (1 − λ)H 0 + λH 1 is the mixed Hessian given in terms of the Hessians H 0 and H 1 , all evaluated at the MECI.The vibrational contributions to the partition functions Z 0 and Z ‡ are the classical harmonicoscillator partition functions of the reactant and MECI, where the latter can be evaluated by computing the frequencies of the mixed Hessian H after projecting out the reaction coordinate with P = I − ∆G ⊗ ∆G T .Here I is the f -dimensional identity matrix and ∆G

d 2 S
dτ 2 = − d 2 S dt 2 and defined S = d 2 S dτ 2 , which is negative at the instanton.As shown in Ref. 15, this derivation can equally be carried out in the ring-polymer formalism and is easily amended to account for rotational and translational contributions.In this ring-polymer version of the instanton rate, presented in Eq. (3.3) of Ref. 15, χ corresponds to the (r 0,c , r N R ,c )-element of the inverse of the second-derivative matrix of the action [defined in Eq. (3.2) of Ref. 15] with respect to all bead positions multiplied by α 2 .Note that in Ref. 15, the reactant and product states were labelled by n ∈ {R, P}, whereas in this work we use n ∈ {0, 1}.The notation of Ref. 15 uses the notation r 0,c , r N R ,c to refer to the coupling coordinate of the hopping beads, x c , with i = 0 and i = N R .
) of the main text], which may lead to the occurrence of several stationary paths of the effective action, S eff = S − ln |∆(x )∆(x )|.This is illustrated in Fig. SI.2, where we plot the integrand of Eq. (SI.16) for the three-mode model defined in Sec.SI.1 in the x c x c coordinate space.From the figure, we can identify two positive lobes, corresponding to the direct paths that bypass the CI and contribute to the total rate with a positive sign, as well as two negative lobes, corresponding to the winding pathways that give rise to a negative contribution to the rate (see Fig. SI.5 and Fig. 1 of the main text for the pathways).Note that the reason for the negative contribution associated with the winding paths is that they employ two different hopping points, x = x , on opposite sides of the CI.

Figure
Figure SI.2:Colour map of the integrand in Eq. (SI.16) for the model system defined in Sec SI.1 in atomic units.The two positive lobes correspond to the direct paths in the SDI formalism, while the two negative lobes correspond to the winding paths (see also Fig. SI.5).To construct the plot, we transformed the system to a spin-boson model, as described in Eq. (SI.6), and used the corresponding analytic expressions for the action.13The two remaining coordinates, as well as the value of τ were fixed to their values at the stationary point of the (effective) action.Note that in this symmetric spin-boson model, the stationary points in the SCI and SDI formalisms only differ in the coupling coordinate, while the values of the two other modes and τ are the same.

Fig. SI. 3 Figure SI. 3 :
Fig. SI.3(a) at the stationary value of τ = β /2.In previous instanton studies on reactions without CIs,19 the GR flux correlation function at the stationary value of τ was shown to be almost perfectly described by a Gaussian.This is not the case in Fig. SI.3(a), which exhibits dips to negative functional values.This is caused by the two-point correlation, χ(τ + it), which fluctuates relatively quickly in time.The SCI approximation of the correlation function treats χ as a constant and thus overestimates the rate.The SDI approach captures the negative contributions to the correlation function by virtue of locating several instantons with positive and negative contributions to the rate, see also Fig. SI.3(b).This leads to the improved agreement with the GR rate observed inTable SI.2.
Fig. SI.3(a), is in excellent agreement with the quantum-mechanical GR correlation function.With knowledge of χ(τ ), we can approximate the exponential in Eq. (SI.18) by a Gaussian in the time domain and take the final integral over time analytically.

Figure SI. 4 :
Figure SI.4:Rates computed with different methods for the three-mode model system.
Fig. SI.4,despite the neglect of tunnelling.The SDI instantons are depicted in Fig. SI.5 along with the stationary points of V eff projected onto the two-dimensional plane spanned by x t and x c .The Figure confirms that the two stationary points of the effective potential correspond to the classical limit of the curving or direct instantons (blue and red paths), whereas the pair of instantons winding around the CI (violet paths, on top of each other, winding clockwise and counterclockwise), which reduce the total rate, do not have a classical analogue.
Fig. SI.6, we compare our first-principles results based on CDFT with rate constants computed for the full-dimensional LVC model constructed in Ref. 25.It can be seen that the SCI result for the LVC model at 300 K is in agreement with the GR rate constant extracted from the short-time dynamics shown in Fig. 18 of Ref. 25.This is further evidence that the semiclassical approximation inherent in instanton theory is well-suited for the description of the charge transfer in the BMA cation.While the first-principles results agree with the rate constants for the LVC model at high temperatures, a significant deviation of more than one order of magnitude is observed in Fig. SI.6 at lower temperatures.One reason for this difference is that different electronic-structure methods were used in the two studies.In particular, CASSCF was used to construct the LVC model, whereas our first-principles calculations are

Figure SI. 6 :
Figure SI.6:Rate constants computed with SCI for the charge transfer in the BMA cation based on on-the-fly CDFT evaluations of the PESs.For comparison, we present results for the full-dimensional LVC model of Ref. 25.The "LVC FGR" results are extracted from Fig. 18 of Ref. 25, and the "LVC SCI" results are computed with instanton theory on the same LVC model.

Table SI . 1 :
Parameters of the three-mode model.

Table SI .
2 seems to suggest that the classical rate is in reasonably good agreement with the exact result.This is simply a fortuitous cancellation of errors, on which one cannot rely in general.

Table SI .
3: Rate constants (in units of ns −1 ) for the charge transfer process in BMA cation (top) and deuterated BMA cation (bottom) computed with CDFT.Table SI.4:Same as Table 2 in the main text but for fully deuterated BMA.

Table SI .
6: Same as Table2in the main text but for the LVC model..A second difference is that the LVC model uses a global harmonic approximation, whereas our first-principles instanton calculations make no assumptions about the shape of the underlying PESs.The comparison between the different contributions to the SCI rates based on our onthe-fly CDFT calculations (Table2of the main text) and the LVC model (TableSI.6)