Solvation or not solvation: tunneling reactions of molecules embedded in cryogenic matrices

Abstract

Noble gas cryogenic matrices have been widely employed to study the tunneling reactivity of embedded molecules. We regard cryogenic matrices as an elementary type of solvent, because the interactions between solvent molecules and between the solvent and the solute are limited to weak, non-directional forces. Nevertheless, they remain very interesting for fundamental solvation studies, as the tunneling reactivity varies considerably for different noble gases. Here, we present a multidimensional nuclear quantum approach based on the semiclassical transition state theory approximation with a quantum mechanics/molecular mechanics potential to investigate the effect of different cryo-matrix environments on the H-tunneling CO bond rotamerization reaction. Specifically, we show how the rotamerization reaction rate of glycine, formic acid, acetic acid, and their deuterated variants changes with different types of cryo-matrix embedding and in the gas phase. After reproducing the available experimental results, we show that cryogenic matrices indeed play a crucial role in the tunneling process, which is not limited to washing out or quenching tunneling. We conclude that condensed phase matrices can enhance the reaction rate by interacting with substituents at the α-carbon site when present. Our approach opens the possibility for future studies of more complex solvation scenarios to gain physical insight into the effect of solvation on tunneling-dominated reactions.

---

1 Introduction

Understanding the effects of solvation on chemical reactivity is a pivotal task.^1–5 Solvents usually have a complex internal structure and interaction networks; therefore, investigations based on a gradual approach in which the solvent is as elementary as possible are ideal for fundamental solvent effect studies. For this reason, in the present work, we focus on the conformational tunneling reactivity of a set of similar organic acids in noble gas cryo-crystals. We regard cryogenic noble gas matrices as a simple solvent prototype,⁶ where the solvent is composed of single atoms interacting through weak and non-directional forces both with each other and with the solute. This situation is simpler than other solvents, such as water, which can establish solvent structures and hydrogen bonds with the solute and among themselves. The relative simplicity of cryogenic noble gas matrices will facilitate modeling and thus understanding the physical picture related to solvent effects on reactivity. Moreover, experimental half-life times and rates are available for the reactions under examination, as well as X-ray diffraction experimental structures, which provide the actual geometry of the solvent.⁷

The goal of the present work is to assess how noble-gas cryogenic matrices are relevant for the acid conformational reactivity, i.e., how these matrices play a role in the reactivity by acting as a solvent. From a general perspective, the solvent changes the energy of the transition state (TS) of the reaction due to its polarity, a bulk effect, or by establishing specific molecular interactions between the TS structure and one or more solvent molecules.⁶ Addressing the influence of the cryogenic matrix environment on reactive processes involving quantum tunneling raises several key questions: How does the environment influence the tunneling mechanism? Is it a bulk or local effect? How do we model the action of the solvent when quantum effects play a role? And how can we reliably calculate tunneling rates by including the environment in a computationally feasible way? The answer to these questions will help to understand how solvents can alter the selectivity of reactions driven by nuclear quantum effects, such as hydrogen or heavy-atom tunneling.^6,8–11 This understanding is particularly important for tunneling mechanisms in biologically relevant molecules.^12–17 The investigation of these scenarios will become manageable with our computationally feasible approach, which accurately incorporates solvent effects into tunneling rate calculations.

Experimentally, cryo-matrix embedding is well-established,^18–27 and is commonly used to approximate the reactivity of isolated molecules. It has often been assumed that gas-phase theoretical results can be directly compared with these experiments. This comparison would be strictly valid in the absence of interactions between different matrix reactive sites and between noble gas atoms and the embedded molecule.^28–31 However, it has been largely demonstrated that even noble gases can significantly interact with simple molecules, impairing the theoretical evaluation of rates using the gas phase approximation.^32–35 For these reasons, there have been developments to model the environment effect on tunneling reactions by accounting for the change in the reaction barrier upon solvation or by hypothesizing empirical corrections due to changes in the polarizability of the environment.^8,9 Within this approach, solvated systems are generally modeled using implicit solvent methods, i.e., without atomistic resolution, or by selecting a single or a few most important solvent molecules involved in the interaction with the reacting moiety.^36–39 Then, one can resort, for example, to the classical transition state theory (TST) with tunneling corrections. Usually, one-dimensional tunneling models are employed, assuming the reaction coordinate to be separable. Additionally, semi-empirical models that utilize the Einstein and Debye approximations have been recently employed for reactions occurring in cryogenic matrices. They can predict the temperature dependence of the rate, but still require experimental data for fitting.⁴⁰ More generally, accounting for the influence of the matrix environment is not straightforward; to obtain experimentally accurate results, a multidimensional nuclear quantum approach may be necessary, as we show in the present work.⁴¹

2 Theory

Semiclassical transition state theory (SCTST) is a quantum TST that accounts for multidimensional tunneling by considering anharmonicity and couplings of the reaction coordinate and of the bound vibrational modes.^42,43 SCTST is ideal for this application, as we can account for the influence of matrix phonon modes on the tunneling reactivity. Pioneered by W. H. Miller,⁴² the SCTST formulation of the rate constant iswhere Q^tra_TS(R)(T) and Q^rot_TS(R)(T) are the translational and rotational partition functions for the transition state (reactant). In our calculations we include anharmonic effects since Q^vib_R(T) is the fully coupled second order vibrational perturbation theory (VPT2) vibrational partition function for the reactant molecule, and it is recovered from the reactant vibrational density of states (DOS) ρ(E) as N^SC(E) is the cumulative reaction probability (CRP) in the semiclassical approximation.^42,44–46 Calculating the SCTST rate constants is computationally cost-effective because it involves static information only, such as the energy barrier that separates reactants from products, the harmonic frequencies of the reactant and the transition state, and the anharmonic couplings determined using VPT2 approximation.⁴⁷

Our group has recently developed and applied SCTST to a wide range of problems, from the investigation of the heavy atom tunneling phenomena up to room temperature in organic reactions^30,48 to the accurate estimate of cryogenic reaction rates for isomer interconversion reactions in the gas phase using high-level ab initio methods such as coupled cluster CCSD(T)-F12b.⁴⁹ To enable these computations and overcome the curse of dimensionality, fast, parallel and efficient SCTST codes for the evaluation of the density of states,⁵⁰ the CRP,⁵¹ and the anharmonic constants⁴⁸ have been developed and this allows us now to extend the applicability of SCTST to solvated systems.

This work focuses on the rotamerization reactions of the OH group from cis to trans configurations in small organic molecules containing the carboxylic acid moiety. H-atom tunneling dominates the isomerization processes at cryogenic temperatures, and the experimental reaction rates show a dependence on the gas matrix employed to embed and isolate the molecules. The collection of reactions that we are simulating is reported in Fig. 1. These are the glycine (GLY),³⁴ the formic acid (FA),^33,52 and the acetic acid (AA) molecules⁵³ OH bond cis to trans isomerization. Considering that we are focusing on tunneling reactions, we included the tri-deuterated glycine (TDGLY) and the methyl-deuterated acetic acid (DAA), where the deuteration sites are reported in yellow in Fig. 1.

Fig. 1 (a) Glycine interconversion reaction. (b) Tri-deuterated glycine reaction. (c) Formic acid interconversion reaction. (d) Acetic acid reaction. (e) Methyl-deuterated acetic acid interconversion reaction. Deuterium is represented in yellow.

To simulate the explicit cryo-matrix embedding of the reactive molecule and the influence of the bulk effect on the reactivity, here we propose a new approach that combines SCTST and our own N-layered integrated molecular orbital and molecular mechanics (ONIOM) protocol.^54,55 Specifically, we employ a 2-layer ONIOM among the possible multiscale methods. This implementation provides a simple and effective way to account for matrix effects in complex systems. The system is divided into two parts. The reactive molecule is treated at the quantum mechanical (QM) level. We chose the 2nd-order Moller–Plesset perturbation theory (MP2) with the aVDZ basis set,⁵⁶ which proved to be sufficiently accurate for estimating the energetics of these reactions.³⁷ We have also performed density functional theory (DFT) calculations using the B3LYP functional with the def2-TZVP basis set with Grimme dispersion treatment D3 and BJ damping.^57–59 However, we found a consistent overestimation of reaction rate constants with DFT (see Fig. S2 and S3), so we relied on the higher-level MP2 calculations. To model the environment surrounding the reactive center at the cheaper molecular mechanics (MM) level, we employed the universal force field (UFF).⁶⁰ The ONIOM-UFF setup has been proven to be successful in treating different types of chemical problems.⁵⁵ The detailed setup steps of the ONIOM-SCTST protocol are described in the SI.

The size of the noble gas lattice, the cavity dimension, and the embedded molecule geometry are determined through geometry optimization. In the FA case, the starting cavities were obtained from the literature.³⁶ For GLY and the AA, the same initial cavities were used as a guess and then modified by removing additional noble gas atoms to fit the new molecules. The right panel of Fig. 2 shows an example of the krypton-embedded glycine (GLY-Kr) molecule.

Fig. 2 (Right panel) Glycine molecule embedded in the cryogenic Kr gas matrix cluster (GLY-Kr) upon geometry optimization at the TS state geometry with the ONIOM(MP2/aug-cc-pVDZ:UFF) level of theory. (Left panel) The two low-frequency normal modes (NM) that are mostly coupled to the modes of the molecule, and their phonon frequency. For each of them, we show in color the nuclei that exhibit significant normal-mode displacements among all the atoms in the structure and the displacement directions (blue vectors). Displacements of module lower than 0.0016 au are not visualized.

In our optimization process, we first found the number of noble-gas atoms sufficient to reproduce a cryo-crystal. The size of the simulated system was determined by keeping the molecule and the cavity fixed and extending the lattice sizes until the optimization no longer perturbed the geometry of the atoms at the edges of the matrix cluster. Then, the cavity is optimized by exploiting the calculation of the formation energy of the system, which accounts for the interaction energy and the sum of the deformation energies for the lattice and the embedded molecule.³⁶ We report the final optimized geometries in the SI. Upon optimization, the full-dimensional Hessian is calculated using the available ONIOM second derivatives,⁵⁵ and after diagonalization, the vibrational normal modes are found. We then identify a set of modes (active modes) by selecting those whose displacements are mainly localized on the QM molecule, in addition to possible other normal modes that involve contributions from the noble gas matrix displacements, while showing a relevant contribution to the central molecule. A couple of the selected active modes involving normal-mode displacements on both the embedded molecule and some of the matrix atoms are reported in the left panel of Fig. 2 for the GLY-Kr system. We report the active mode list for all systems in Tables S2–S4 of the SI. We perform anharmonic VPT2 calculations on the active modes only. However, the energies used to compute the high-order derivatives are full-dimensional, including the contribution from all atoms in the system. For the ONIOM calculations, we employed the Gaussian16 code,⁶¹ and we implemented the necessary features to pair it with the FDACC code that takes care of the parallel VPT2 calculation.^48,62

3 Results

We start from reaction (a) of Fig. 1 and report in Fig. 3 our ONIOM-SCTST GLY interconversion half-life times, where the QM is at the level of Moller–Plesset (MP2) theory. For the gas phase, we also report the CCSD(T)-F12b results. The available experimental values at 12 K are shown for comparison. In each plot, the gas-phase results are reported in violet. The MP2 gas-phase half-life time estimates are consistently close to the CCSD(T)-F12b ones, indicating that MP2 can accurately describe the system, and we can safely employ it for more computationally intensive condensed-phase calculations. For each matrix, either Ar, Kr, or Xe, the half-life time is significantly reduced over the entire range of cryogenic temperatures compared to the gas phase, suggesting a significant solvation effect of the matrices. In particular, the calculations quantitatively reproduce the experimentally observed half-life time at 12 K for all matrices.

Fig. 3 Half-life time calculations for the glycine interconversion in noble gas matrices (reaction (a) of Fig. 1). Each panel is for a different matrix. Panel (a) for argon cryo-matrix, (b) for krypton and (c) for xenon ones. Vertical bars are the experimental estimates with their error bar.³⁴ SC/F12-CCSD(T) is the SCTST calculation at the coupled cluster level of theory, while SC/MP2 is at the Moller–Plesset level.⁴⁹

When the tunneling H-atom and other hydrogen atoms are deuterated, as in the TDGLY reported in Fig. 1(b), the half-life time is orders of magnitude longer. Also in the deuterated case, the reduction of the half-life time in the matrix with respect to the gas-phase one is on average a factor of 3, which is comparable with the hydrogenated cases. Once again, the good agreement among different post-HF methods is confirmed in the gas phase calculations reported in Fig. 4. Specifically, the CCSD(T) and MP2 results are in agreement with the CCSD(T)-F12b ones. MP2 results are again a good estimate of the experimental values in the condensed phase.

Fig. 4 Tri-deuterated glycine interconversion half-life times in noble gas matrices (reaction (b) on Fig. 1). The black symbols refer to the experimental values at 12 K,³⁴ the continuous lines represent the SCTST gas phase results at the coupled cluster level, and the dashed lines the SCTST ones at the MP2 ab initio level of theory.

A completely different picture is observed for the FA. Even if the FA molecular structure closely resembles that of GLY, and the interconversion process involves the same OH moiety and matrix environment, we observe different half-life times (and therefore different reaction rate constants as reported in Fig. S1 in the SI) for different matrices, as already shown by experiments.⁶³

Our SCTST results for the FA interconversion are compared directly with the experimental ones³³ in Fig. 5, where the original experimental reaction rates k(T) have been converted into half-life times t_1/2 using the relation t_1/2 = ln 2/k(T). For the reaction in Xe, the experiments show some discrepancies, and two different sets of values are reported in the literature for the same temperature range. Overall, the SCTST results are accurate except for the Xe matrix case in the very low temperature regime. When comparing the FA interconversion with the GLY one, as shown in Fig. 5, the FA half-life times are longer in the condensed phase rather than in the gas phase, and this is at odds with what we have found for GLY and TDGLY.

Fig. 5 Formic acid interconversion half-life times in noble gas matrices (reaction (c) in Fig. 1). The diamonds represent the experimental results.³³ The violet dashed lines in the three panels refer to the SCTST gas-phase calculations.

To gain a deeper insight into this opposite trend, we consider that even if the reactions presented here involve the same isomerization of the carboxylic acid moiety, the reaction coordinate for the GLY molecule interconversion involves hydrogen atoms on the amino group, in addition to the rotation of the OH group along the CO bond axis. Indeed, the displacements of the reactive normal mode indicate a twisting motion of hydrogen atoms of the amine moiety during the reaction, as one can appreciate from panel (b) of Fig. S10. To understand the role played by the presence of an α-carbon and its substituents in terms of reactivity, we analyzed the vibrational frequencies for each noble-gas matrix. Table S2 presents the normal mode frequency values for the GLY reaction in different cryogenic matrices for the reactant and transition state geometries. For reference, we reported the same data in Table S3, but for the FA reaction. A comparison between the two tables shows that FA is characterized by much higher frequency values than GLY. Specifically, the percentage changes that the vibrational frequencies undergo for each environment are definitely smaller for the FA than for GLY. In other words, the interconversion in the FA molecule is less affected by the presence of the cryogenic matrix. This is an atomistic perspective on the different reactive tunneling behaviors observed above, and it suggests, as part of the explanation, that the primary solvent effect is the bond stiffening predominantly caused by electron cloud repulsions, and the consequential changes in the vibrational frequency values. Importantly, this phenomenon appears to be sensitive to the type of substituents on the α-carbon. Deviations from gas-phase vibrational frequencies become more pronounced for α-substituted molecules, proving that this part of the molecule is interacting with the matrix. This explains why GLY and TDGLY behave differently by changing the solvent, and that the FA, which shows low-frequency modes at much higher values than GLY and TDGLY, i.e., at about 600 cm⁻¹, exhibits reduced interaction with the surrounding matrix. We deduce that when the molecule presents low-frequency modes comparable to those of the matrix, these modes mix with the matrix and act essentially as conduits, facilitating energy transfer via vibrational coupling to other modes. FA lacks these bridging modes, leading to less efficient vibrational energy transfer.

Along this line, we investigate the half-life time changes when the three methyl hydrogen atoms in the AA molecule are replaced with deuterium, i.e., the DAA reaction of panel (d) Fig. 2.⁵³ This is an example of a secondary isotope effect, as the deuteration involves hydrogen atoms that do not participate directly in the conformational change of the OH group. The atomic displacements involved in the reactive normal mode are presented in Fig. S10 panel (a), and the isomerization involves rotation of the methyl moiety. As reported in Fig. 6, we first reproduce the experimental value for the non-deuterated AA molecule in the Xe matrix with our ONIOM-SCTST MP2 method. Then, we employ this setup to show the AA condensed-phase rate enhancement upon methyl deuteration.

Fig. 6 SCTST half-life times at the MP2/aVDZ level of theory for the interconversion of the acetic acid molecule and its methyl-deuterated isotopologue in xenon cryogenic matrix and gas phase (reactions (d) and (e) of Fig. 1). Black diamond represent the experimental values in xenon.⁵³

Remarkably, we observe that the effect of deuteration is the opposite for the same reaction in vacuum conditions than in matrices, i.e., the reaction in the gas-phase is slower after deuteration, as reported in violet in Fig. 6. In contrast, in the matrix environment, the reaction is an example of an inverse kinetic isotope effect. These results confirm that the presence of the α-carbon and other additional substituents may effectively alter the reactivity in matrices. Since deuteration lowers the frequency values, we hypothesize that this alteration can lead to a more efficient vibrational coupling with the low-frequency phonon modes of the matrix and account for the changes in the rates.

4 Discussion

In light of these results, we can make some considerations about the solvation process for reactions dominated by tunneling. Even if tunneling is an instantaneous process and the shape of the minimum energy path does not play a direct role, the rates of the cryo-matrix-embedded molecule reactions shown above are always different from the gas phase. The solvation effects of the cryogenic matrices accelerate the reaction when substituents are present on the α-carbon (see the GLY, TDGLY, AA, and DAA cases). Otherwise (FA case), the presence of the matrix slows down the interconversion of the OH group, as one would have expected. Additionally, we also observed that for FA, the heavier matrices show stronger interactions. By comparing the GLY and FA reactivities, we deduce that the low-frequency phonon modes of the matrix mix with those of glycine and accelerate the reaction rate.

One may wonder if we really need such atomistic detail to achieve this level of accuracy. To clarify this question, we conducted further investigations towards a better understanding of the effects of electronic cloud repulsions through implicit solvation models. This investigation emphasizes the importance of dispersion interaction modeling and the significant role of explicitly incorporating the noble gas matrix when simulating chemical kinetics. We employed two different implicit solvent approximations. One is the polarizable continuum model (PCM)⁶⁴ and the other is the Universal Solvation Model based on Solute Electron Density (SMD).⁶⁵ Details on the (PCM/SMD)-SCTST anharmonic analysis and rate constants are reported in the SI (see Fig. S4–S9). Specifically, PCM-SCTST falls significantly short in both quantitative and qualitative aspects when predicting reaction rate trends in these solvation environments. In contrast, the SMD model can accurately capture the behaviors of FA and GLY qualitatively across matrices, although we do not achieve quantitative agreement comparable to that of atomistic simulations. We expected that PCM would underperform in noble gas matrices where the predominant interactions with reactive molecules are short-range and non-electrostatic, including electronic exchange repulsion and dispersion forces. Instead, in approaches like the ONIOM-SCTST employing UFF, these interactions are considered, even if generally simplified into a single van der Waals term. As far as the SMD modeling is concerned, the SMD-SCTST methodology effectively accounts for short-range forces, aligning well with the experimental qualitative observations. Nevertheless, the SMD-SCTST rates are much less accurate than the ONIOM-SCTST in predicting the experimental rates (see SI Fig. S4–S9), suggesting that the directionality of the charge fluctuation and the atomistic detail in matrix interactions are pivotal.

One can also try to rationalize our results in terms of barrier height and frequency in a one-dimensional picture along the reactive path. Table S1 of the SI reports these quantities for the GLY, FA, and AA systems for all the cryo reactions reported above. In Table S1, in addition to the classical barrier, we consider under the section “Anharmonic Barrier” also the barrier height corrected by the difference in the anharmonic Zero Point Energy (ZPE) of the reactants and the transition state. By inspection of either the classical barrier or the anharmonic barrier values, one would deduce that GLY OH rotamerization would be faster in the cryo-matrix environments than in the gas phase, and the opposite for the FA. This is in agreement with our findings. However, one should consider that the OH rotamerization is essentially a tunneling reaction in this case, and the barrier width plays a very important role, as extensively discussed in ref. 66. The wider the barrier, the smaller the frequency along the reaction path at the transition geometry, and the smaller the tunneling transmission probability. Table S1 reports a trend opposite to what one would expect. This shows that the enhanced rate is not just a matter of tunneling. For both GLY and FA interconversion, Table S1 shows that the barrier width increases when moving from the gas to the condensed phase, indicating a smaller amount of tunneling. Eventually, one would not have been able to draw definitive conclusions based just on the information reported in Table S1, and we confirm once again that the multidimensional anharmonic SCTST calculations are necessary.

Then, one may wonder if it is really necessary to employ a high level of electronic structure theory^67,68 or if the cheaper DFT approach is sufficient. Fig. S2 and S3 address this issue for the GLY and FA reactions, respectively. Fig. S2 clearly shows that the DFT gas-phase setup better reproduces the experimental results than the condensed-phase setup. In the case of the FA reactions, reported in Fig. S3, the agreement of the simulations with the cryo-matrix results is indeed improved, but almost an order of magnitude less accurate than the MP2 ones.

5 Conclusions

In conclusion, we show that a high ab initio level of theory and a multidimensional SCTST tunneling rate calculation allow us to obtain reaction rate constants accurate enough to appreciate the role of the cryogenic matrices. Inclusion of nuclear quantum effects, anharmonicities, and atomistic solvent details in the computation have been necessary to achieve the goal of reproducing the experiments. This has also been shown recently for IR simulations of relevant biological molecules in water.^69,70 Our findings indicate that short-range exchange repulsion and dispersion forces influence the reactivity of molecules within these cryogenic environments, with the extent of this influence depending on the specific noble gas matrix employed. These conclusions are in agreement with previous spectroscopic investigations in similar systems.²⁰ The solvent or matrix plays a quite active role as confirmed, and it can either enhance or decrease the reaction rate.⁶ Thus, understanding the solvation process for tunneling reactions is not trivial even at low temperatures in cryogenic matrices, as we observe opposite trends for the same type of interconversion and for different noble-gas matrices. If one thinks that the solvent degrees of freedom might only dissipate the reactant energy, providing the reactants with additional routes and, therefore, decreasing the rate of these tunneling transmission, in the present work, we have shown that this is not the case. The comparison between gas-phase and condensed-phase rates clearly suggests that the solvent plays the active role of enhancing the rate. This is a case where the solvent does not act as a friction source that suppresses quantum tunneling in a condensed phase environment.⁷¹ Specifically for the H-tunneling CO bond rotamerization reaction, we find that substituents at the α-carbon mostly interact with the low-frequency modes of the bulk, accelerating the reaction rate.

Our ONIOM-SCTST approach offers a practical method for accounting for the effects of cryogenic gas matrices on the reactivity of embedded molecules, successfully narrowing the gap between experimental results and theoretical predictions for such reactions. The ONIOM-SCTST method paves the way for exploration of complex solvation environments. These include molecules such as nitrogen, which is known to interact strongly with the carboxylic acid group, and potentially other solvation scenarios with polar molecules.^40,72 Addressing these situations would require different force fields to accurately account for polarizability and intramolecular vibrations, such as AMOEBA FF.⁷³

Author contributions

G. Mandelli: data curation (supporting), formal analysis (lead), investigation (lead), methodology (equal), software (lead), validation (equal), visualization (equal), writing – original draft (equal), writing – review & editing (supporting). C. Aieta: conceptualization (lead), data curation (supporting), formal analysis (supporting), investigation (supporting), methodology (equal), project administration (equal), software (supporting), supervision (equal), validation (equal), visualization (equal), writing – original draft (equal), writing – review & editing (lead). M. Ceotto: conceptualization (supporting), data curation (lead), funding acquisition (lead), methodology (supporting), project administration (equal), resources (lead), supervision (equal), validation (equal), writing – original draft (equal), writing – review & editing (supporting).

Conflicts of interest

The authors declare no conflict of interest.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5sc06818e.

Acknowledgements

The authors thank Prof. John Barker, Prof. Riccardo Conte, Dr Eno Paenurk, and Luca Corneo for the useful discussions. G. M. thanks the Università degli Studi di Milano for a PhD scholarship. Part of the needed CPU time was provided by CINECA (Italian Supercomputing Centre) under the ISCRAC project SCTSTGly – HP10C0CQE0.

Notes and references

1. C. Reichardt, Org. Process Res. Dev., 2007, 11, 105–113.
2. A. Gutberlet, G. Schwaab, Ö. Birer, M. Masia, A. Kaczmarek, H. Forbert, M. Havenith and D. Marx, Science, 2009, 324, 1545–1548.
3. R. Jimenez, G. R. Fleming, P. Kumar and M. Maroncelli, Nature, 1994, 369, 471–473.
4. V. Conti Nibali and M. Havenith, J. Am. Chem. Soc., 2014, 136, 12800–12807.
5. B. Mennucci, J. Phys. Chem. Lett., 2010, 1, 1666–1674.
6. T. Schleif, M. Prado Merini, S. Henkel and W. Sander, Acc. Chem. Res., 2022, 55, 2180–2190.
7. R. Wyckoff, in Crystal Structures, Interscience Publishers, New York, London, Sydney, vol. 1, 1963.
8. T. Schleif, M. Prado Merini and W. Sander, Angew. Chem., Int. Ed., 2020, 59, 20318–20322.
9. T. Schleif, J. Tatchen, J. F. Rowen, F. Beyer, E. Sanchez-Garcia and W. Sander, Chem. – Eur. J., 2020, 26, 10452–10458.
10. T. Schleif, J. Mieres-Perez, S. Henkel, M. Ertelt, W. T. Borden and W. Sander, Angew. Chem., Int. Ed., 2017, 56, 10746–10749.
11. Y. Zhou, W. Fang, L. Wang, X. Zeng, D. H. Zhang and M. Zhou, J. Am. Chem. Soc., 2023, 145, 8817–8821.
12. J. P. Layfield and S. Hammes-Schiffer, Chem. Rev., 2014, 114, 3466–3494.
13. T. Ichikawa, T. Takahashi, H. Koizumi and T. Takada, J. Phys. Chem. A, 2000, 104, 2581–2586.
14. P. R. Schreiner, Trends Chem., 2020, 2, 980–989.
15. B. A. Johnson, Y. Hu, K. Houk and M. A. Garcia-Garibay, J. Am. Chem. Soc., 2001, 123, 6941–6942.
16. E. Pollak, J. Phys. Chem. B, 2012, 116, 12966–12971.
17. J. P. Roque, C. M. Nunes, E. Ntungwe, A. Fernández-Ramos and R. Fausto, J. Org. Chem., 2025, 90, 10476–10486.
18. R. T. Hall and G. C. Pimentel, J. Chem. Phys., 1963, 38, 1889–1897.
19. M. Pettersson, J. Lundell, L. Khriachtchev and M. Räsänen, J. Am. Chem. Soc., 1997, 119, 11715–11716.
20. R. Fausto, G. O. Ildiz and C. M. Nunes, Chem. Soc. Rev., 2022, 51, 2853–2872.
21. C. M. Nunes, J. P. Roque, S. Doddipatla, S. A. Wood, R. J. McMahon and R. Fausto, J. Am. Chem. Soc., 2022, 144, 20866–20874.
22. A. J. Lopes Jesus, C. M. Nunes and I. Reva, Photochem, 2022, 2, 405–422.
23. J. P. Roque, C. M. Nunes, L. P. Viegas, N. A. Pereira, T. M. Pinho e Melo, P. R. Schreiner and R. Fausto, J. Am. Chem. Soc., 2021, 143, 8266–8271.
24. C. M. Nunes, L. P. Viegas, S. A. Wood, J. P. Roque, R. J. McMahon and R. Fausto, Angew. Chem., Int. Ed., 2020, 59, 17622–17627.
25. G. Bazsó, E. E. Najbauer, G. Magyarfalvi and G. Tarczay, J. Phys. Chem. A, 2013, 117, 1952–1962.
26. E. E. Najbauer, G. Bazso, S. Gobi, G. Magyarfalvi and G. Tarczay, J. Phys. Chem. B, 2014, 118, 2093–2103.
27. V. D. Drabkin and A. K. Eckhardt, J. Phys. Chem. Lett., 2025, 16, 2223–2230.
28. Z. Wu, R. Feng, H. Li, J. Xu, G. Deng, M. Abe, D. Bégué, K. Liu and X. Zeng, Angew. Chem., 2017, 129, 15878–15882.
29. O. Kirshenboim, A. Frenklah and S. Kozuch, Chem. Sci., 2021, 12, 3179–3187.
30. F. Angiolari, G. Mandelli, S. Huppert, C. Aieta and R. Spezia, Chem. – Eur. J., 2024, 30, e202401000.
31. F. Angiolari, S. Huppert and R. Spezia, Phys. Chem. Chem. Phys., 2022, 24, 29357–29370.
32. E. M. Maçôas, L. Khriachtchev, M. Pettersson, R. Fausto and M. Räsänen, J. Phys. Chem. A, 2005, 109, 3617–3625.
33. M. Pettersson, E. Maçôas, L. Khriachtchev, J. Lundell, R. Fausto and M. Räsänen, J. Chem. Phys., 2002, 117, 9095–9098.
34. G. Bazso, G. Magyarfalvi and G. Tarczay, J. Phys. Chem. A, 2012, 116, 10539–10547.
35. I. Sedgi and S. Kozuch, Phys. Chem. Chem. Phys., 2020, 22, 17725–17730.
36. S. Stepanian and L. Adamowicz, Low Temp. Phys., 2020, 46, 155–164.
37. L. I. Trakhtenberg, A. A. Fokeyev, A. S. Zyubin, A. M. Mebel and S. Lin, J. Phys. Chem. B, 2010, 114, 17102–17112.
38. J. Dokic, M. Gothe, J. Wirth, M. V. Peters, J. Schwarz, S. Hecht and P. Saalfrank, J. Phys. Chem. A, 2009, 113, 6763–6773.
39. A. K. Eckhardt, D. Gerbig and P. R. Schreiner, J. Phys. Chem. A, 2018, 122, 1488–1495.
40. K. Marushkevich, M. Rasanen and L. Khriachtchev, J. Phys. Chem. A, 2010, 114, 10584–10589.
41. W. H. Miller, Acc. Chem. Res., 1976, 9, 306–312.
42. R. Hernandez and W. H. Miller, Chem. Phys. Lett., 1993, 214, 129–136.
43. T. L. Nguyen, J. R. Barker, and J. F. Stanton, Atmospheric reaction rate constants and kinetic isotope effects computed using the HEAT protocol and semi-classical transition state theory, in Advances in Atmospheric Chemistry, ed. J. R. Barker, A. L. Steiner, and T. J. Wallington, World Scientific, 2017, ch. 6, vol. 1, pp. 403–492.
44. W. H, R. Hernandez, N. C. Handy, D. Jayatilaka and A. Willetts, Chem. Phys. Lett., 1990, 172, 62–68.
45. W. H. Miller, J. Chem. Phys., 1974, 61, 1823–1834.
46. W. H. Miller, J. Chem. Phys., 1975, 62, 1899–1906.
47. A. M. Rosnik and W. F. Polik, Mol. Phys., 2014, 112, 261–300.
48. G. Mandelli, C. Aieta and M. Ceotto, J. Chem. Theory Comput., 2022, 18, 623–637.
49. G. Mandelli, L. Corneo and C. Aieta, J. Phys. Chem. Lett., 2023, 14, 9996–10002.
50. C. Aieta, F. Gabas and M. Ceotto, J. Phys. Chem. A, 2016, 120, 4853–4862.
51. C. Aieta, F. Gabas and M. Ceotto, J. Chem. Theory Comput., 2019, 15, 2142–2153.
52. L. I. Trakhtenberg, A. A. Fokeyev and A. M. Mebel, Chem. Phys. Lett., 2013, 574, 47–50.
53. E. Maçôas, L. Khriachtchev, M. Pettersson, R. Fausto and M. Räsänen, J. Chem. Phys., 2004, 121, 1331–1338.
54. S. Dapprich, I. Komáromi, K. S. Byun, K. Morokuma and M. J. Frisch, J. Mol. Struct.:THEOCHEM, 1999, 461, 1–21.
55. L. W. Chung, W. Sameera, R. Ramozzi, A. J. Page, M. Hatanaka, G. P. Petrova, T. V. Harris, X. Li, Z. Ke and F. Liu, et al., Chem. Rev., 2015, 115, 5678–5796.
56. R. A. Kendall, T. H. Dunning Jr and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
57. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
58. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
59. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
60. A. K. Rappé, C. J. Casewit, K. Colwell, W. A. Goddard III and W. M. Skiff, J. Am. Chem. Soc., 1992, 114, 10024–10035.
61. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian 16 Revision C.01, Gaussian Inc., Wallingford CT, 2016.
62. J. R. Barker, MultiWell-2023 Software Suite, University of Michigan, Ann Arbor, Michigan, USA, 2023, http://clasp-research.engin.umich.edu/multiwell/.
63. J. P. Wagner, H. P. Reisenauer, V. Hirvonen, C.-H. Wu, J. L. Tyberg, W. D. Allen and P. R. Schreiner, Chem. Commun., 2016, 52, 7858–7861.
64. B. Mennucci, M. Cossi and J. Tomasi, J. Chem. Phys., 1995, 102, 6837–6845.
65. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396.
66. G. Qiu and P. R. Schreiner, ACS Cent. Sci., 2023, 9, 2129–2137.
67. G. Knizia, T. B. Adler and H.-J. Werner, J. Chem. Phys., 2009, 130, 054104.
68. T. B. Adler, G. Knizia and H.-J. Werner, J. Chem. Phys., 2007, 127, 221106.
69. D. Moscato, G. Mandelli, M. Bondanza, F. Lipparini, R. Conte, B. Mennucci and M. Ceotto, J. Am. Chem. Soc., 2024, 146, 8179–8188.
70. G. Schwaab, R. Pérez de Tudela, D. Mani, N. Pal, T. K. Roy, F. Gabas, R. Conte, L. Durán Caballero, M. Ceotto and D. Marx, et al., Phys. Rev. Lett., 2022, 128, 033001.
71. S. Mousavi and S. Miret-Artés, Ann. Phys., 2018, 393, 76–92.
72. S. Lopes, A. V. Domanskaya, R. Fausto, M. Räsänen and L. Khriachtchev, J. Chem. Phys., 2010, 133, 144507.
73. J. W. Ponder, C. Wu, P. Ren, V. S. Pande, J. D. Chodera, M. J. Schnieders, I. Haque, D. L. Mobley, D. S. Lambrecht, R. A. DiStasio, M. Head-Gordon, G. N. I. Clark, M. E. Johnson and T. Head-Gordon, J. Phys. Chem. B, 2010, 114, 2549–2564.

---