Time-Resolved Chemical Bonding Structure Evolution by Direct-Dynamics Chemical Simulations

Time-Resolved Chemical Bonding Structure Evolution by DirectDynamics Chemical Simulations Mario Piris,*Xabier Lopez, and Jesus M. Ugalde Cite This: J. Phys. Chem. Lett. 2024, 15, 12138−12143 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: Direct-dynamics simulations monitor atomic nuclei trajectories during chemical reactions, where chemical bonds are broken and new ones are formed. While they provide valuable information about the ongoing nuclear dynamics, the evolution of the chemical bonds is customarily overlooked, thus, hindering key information about the reaction mechanism. Here we examine the evolution of the chemical bonds for the three main mechanisms of the F−+ CH3CH2Cl reaction using quasi-classical trajectories for the nuclei, and global natural orbitals for the electrons. Key findings include (i) bimolecular nucleophilic substitution (SN2) resembles a one-step bond breaking and formation process; (ii) the elimination mechanisms (synand anti-E2) feature a sequential two-step process of proton abstraction and Cl−elimination; and (iii) the anti-E2 mechanism is slower, exhibits rebound effects, and gets activated by specific vibrational modes. This study highlights the importance of correctly describing and thoroughly analyzing the dynamical evolution of chemical bonds for chemical reaction mechanistic studies. Chemical reactions constitute (complex) rate processes that involve the concerted time evolution of both the nuclei and the electrons of the molecules engaged in the reaction from reactants to products. From a quantum chemistry perspective, the adiabatic time evolution of the nuclei is customarily carried out in the classical domain, and the electrons are assumed to respond quickly enough to the nuclei’s evolution as to justify taking them under the auspices of the Born−Oppenheimer (BO) approximation. 1 Nonadiabatic time evolution dynamics, which requires the consideration of transitions between different electronic states, can be conveniently handled by potential surface hopping methods. 2 Consequently, the study of chemical rate processes consists of (i) selecting the appropriate quantum mechanical method to treat the electrons; (ii) solving the corresponding Schrodinger equation for the energy and estimating its gradient at the position of each nucleus; (iii) making judicious choices of the initial conditions for the nuclei in the ensemble of trajectories to be propagated in time; (iv) numerically integrating the classical equations-of-motion to determine the time evolution of each nucleus in each calculated trajectory; and (v) transforming the trajectories’ final nuclei coordinates and momenta into properties that may be compared with experimental results. As indicated elsewhere, 3 the properties at stake include bond lengths and angles, vibrational, rotational, and translation energies, quantum numbers for the vibrational and rotational degrees of freedom, the amount of energy in individual molecular degrees of freedom, and scattering angles. These are all properties that naturally derive from the positions and momenta of the nuclei but omit further considerations of the chemical bonding structure evolution along the trajectory. The latter is simply not addressed by most chemical dynamics simulation studies. However, a chemical reaction is conceived as a process in which chemical components are transformed to other chemical components. Since chemists identify components by their chemical bonds, i.e., their electronic structure, the evolution of these bonds during the reaction carries significant chemical information encoded in the concept of the reaction mechanism. This refers to the time-ordered sequence of forming and/or breaking chemical bonds, along with the characterization of all transient chemical components throughout the trajectory from reactants to products. 4 Electronic structure calculations enter into chemical dynamics simulations in steps (i) and (ii) mentioned above. However, these are demanding calculations aimed at not only dealing with static considerations of structure but also Received: October 17, 2024 Revised: November 21, 2024 Accepted: November 27, 2024 Published: November 28, 2024 Letterpubs.acs.org/JPCL © 2024 The Authors. Published by American Chemical Society 12138 https://doi.org/10.1021/acs.jpclett.4c03010 J. Phys. Chem. Lett. 2024, 15, 12138−12143 This article is licensed under CC-BY 4.0 Downloaded via UNIV DEL PAIS VASCO on December 11, 2024 at 14:27:01 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. confronting the dynamics of the replacements of some chemical bonds into newly formed ones that takes place in the course of chemical reactions. A word of caution about their appropriateness is in order for most such calculations do not produce chemically acceptable descriptions of the whole evolution of the making and/or breaking of chemical bonds. 5 Accordingly, we have employed the recently introduced Global Natural Orbital Functional (GNOF) 6 to evaluate the BO energies, as previous studies 7,8 have demonstrated that GNOF yields a favorable balance between static and dynamic electron correlations, resulting in chemically informative descriptions and accurate total energies throughout the entire range of internuclear coordinates, from stable structures to their dissociated asymptotic limits. This is a crucial feature for the present investigation, which focuses on the analysis of chemical bonds as they form and break. Additional information on GNOF and its analytic gradients can be found in Appendix. This paper aims to communicate the feasibility of quasiclassical trajectory (QCT) direct-dynamics simulations, supported by GNOF calculations of the electronic structure, to determine the time evolution of the making and breaking of chemical bonds in reactive processes, which enriches the already available information about nuclear dynamics. 9 We have selected for such a purpose the reaction between the fluoride anion and ethyl chloride, an extensively studied canonical case. 10−18 Theoretical simulations conducted to date, along with their comparison to experimental data, have successfully addressed energyand angle-differential reactive scattering cross sections and have effectively distinguished between the competing nucleophilic substitution (SN2) and base-induced elimination (synand anti-E2) reactions. Additionally, it has been established that the reaction of F−+ CH3CH2Cl is predominantly governed by direct dynamics, except at low collision energies. The anti-E2 mechanism has emerged as the dominant mechanism, although the SN2 and syn-E2 mechanisms gain relevance as the collision energy increases. It is worth noting that the present study aims not at replicating already explored aspects of the title reaction but at providing a better understanding of the operating mechanisms. This will be achieved by a thorough analysis of the time evolution of the chemically active bonds along selected trajectories. Thus, based on the previous results, 10−18 we will choose impact parameters and collision energies such that we can monitor reactive trajectories and observe the time evolution of the chemically active natural orbitals and, consequently, the evolution of the chemical bonding structure along the selected reactive trajectories. Earlier high-level electronic structure calculations 19 have identified three primary reaction channels: (i) bimolecular base-induced elimination, which proceeds via both the anti-E2 and syn-E2 pathways, resulting in the formation of CH2�CH2 + HF + Cl−, and (ii) bimolecular nucleophilic substitution (SN2), where the incoming F−attacks the Cαcarbon atom from either the back or front side relative to Cl, producing CH3CH2F + Cl−. However, the front-side attack has a too-high energy barrier of ∼1.3 eV. Accordingly, we adjusted the impact parameter and the direction of the initial velocity vector of the F−anion to favor the three primary mechanisms: (back-side) SN2, anti-E2, and syn-E2. All calculations were performed using the molecular dynamics module integrated into the DoNOF computer program, 20 which employs Beeman’s algorithm 21 to numerically integrate the classical Newton equations-of-motion for the nuclear coordinates, with forces determined on-the-fly as the gradient of the total BO energy at each nuclear position. The initial conditions for Beeman’s numerical integration are selected to represent the quantum mechanical vibrational and rotational energy levels of the reactants in their ground state. Consequently, our calculations should be classified as quasiclassical trajectory (QCT) BO direct-dynamics simulations. 3 The equilibrium structure of ethyl chloride was determined via a full optimization using the GNOF with the cc-pVDZ basis set. 22,23 The molecule was positioned at the center of the laboratory coordinate system, and the initial nuclear momenta were adjusted according to the equilibrium molecular normal vibrational modes. Previous simulations have identified multiple direct and indirect mechanisms of SN2/E2 reactions. It has been established that both reactions transition from a highly indirect character to a direct one as collision energy increases. Schematic potential energy diagrams for the three relevant channels depict similar double-well potential energy surfaces (PESs), with a transition state connecting the reactant and product complexes. 15 The corresponding submerged barriers are extremely close, highlighting the pivotal role of dynamics in the competition between anti-E2, syn-E2, and SN2 pathways. Consequently, we selected a kinetic energy of Ek= 0.2 eV in the laboratory coordinate system for F−, which corresponds to a low collision energy of 0.154 eV in the center-of-mass coordinate system, where indirect mechanisms dominate. The initial separation between the fluoride anion and the center-of-mass of ethyl chloride was set at 6 Å, and each trajectory was integrated until the separation between the final fragments exceeded 6 Å. A time step of 0.1 fs was utilized. The potential energy profiles for the three studied reactive trajectories with F−translational energy Ek= 0.2 eV and CH3CH2Cl in the ground state can be seen in Figure 1, while a movie of these trajectories is available in the Supporting Information. Based on these data, we can conclude that the reactive collisions occur roughly within the time range of 130− 270 fs. The fast fluctuations of the total potential energy for the three mentioned pathways, shown in Figure 1 at both short and long times, reflect the energy oscillations of the quasiFigure 1. Potential energy profiles of reactive trajectories relative to the energy of the reactants, for F−with a translational energy of 0.2 eV and CH3CH2Cl in the ground state. The time step used was 0.1 fs. The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter https://doi.org/10.1021/acs.jpclett.4c03010 J. Phys. Chem. Lett. 2024, 15, 12138−12143 12139 harmonic vibrational modes of the reactants and products, respectively. Observe that, at intermediate times, the oscillation pattern markedly differs from quasi-harmonic, indicative of chemical bonding rearrangement being taking place. The precise nature of the occurring rearrangement for each pathway will be discussed below by inspecting the time evolution of the relevant bonding natural orbitals. The chemical-bonding structure evolution is reflected in a handful of valence natural orbitals. These orbitals, hereafter termed the “chemically active orbitals”, will be discussed for each of three main reaction mechanisms. Full movies of both the nuclear dynamics and the dynamical evolution of the chemically active orbitals may be found in the Supporting Information video files. SN2 Mechanism. Figure 2 shows the time evolution of the chemically active natural orbitals of the SN2 mechanism of the F−+ CH3CH2Cl →Cl−+ CH3CH2F reaction. Notice that these two orbitals describe the breaking of the Cl−Cαbond and the formation of the Cα−F bond. Inspection of their time evolution reveals that these two processes, Cl−Cα bond breaking and Cα−F bond formation, occur synchronously in a single step, which is completed in ∼100 fs time. Ocular inspection of the frames shown in Figure 2 suggests that the transition state has a lifetime of ∼8 fs. Notice that both the σ(Cl−Cα) and σ(Cα−F) bonds begin their breaking and forming processes, respectively, at t= 190 fs and finish at t= 198 fs. syn-E2 Mechanism. Figure 3 shows the orbital time evolution of the chemically active natural orbitals of the syn-E2 mechanism for the F−+ CH3CH2Cl →FH + C2H4+ Cl− reaction. Here we can see the initial proton abstraction from Cβby the incident F−anion as reflected in the time evolution of the two natural orbitals shown in the bottom rows of Figure 3, between t= 160 and 194 fs. Thus, one can clearly observe the formation of the σ(F−H) bond (as seen in the bottom row frames) and the emergence of a carbanion at Cβ, illustrated in the middle row frames. The end result is the formation of a transient hydrogen fluorine “stabilized” chlorinated carbanion depicted in the frame corresponding to t= 194 fs, where these three salient elements of its electronic structure are neatly shown. Namely, the top frame shows the σ(Cl−Cα) bond, the middle frame shows the ClH2Cα−Cβ ⊖H2carbanion, and the bottom frame shows the σ(H−F) bond. From this “critical point” onward, the evolution of the electron density, i.e., the chemical bonding structure, moves on from the vicinity of Cβto Cα. Indeed, inspection of frames t= 200 fs −237 fs reveals the synchronous formation of the new π(C−C) bond and the detachment of the Cl−anion upon breaking of the σ(Cl−Cα) bond. This process takes ∼43 fs. The main observation here is the sequential occurrence of two processes. The former constitutes the proton abstraction by F−from the H−Cβσ-bond and the generation of a carbanion at Cβby retention of the bonding electron pair of the σ(H−Cβ) bond at the Cβcarbon atom. The latter process consists of the rearrangement of the bonding structure around the Cαcarbon to yield a newly formed π(C−C) bond and a chlorine anion which takes the electron pair of the σ(Cl−Cα) bond as it breaks apart. These two processes are decoupled in both timeand realspace coordinate domains. They occur at different reaction times and in different locations. All in all, the above findings Figure 2. Time evolution of the two chemically active natural orbitals for the SN2 mechanism. The bottom line indicates the time stamps (femtoseconds) of the selected snapshots. Contour value = 0.1 au. The top row natural-orbitals’ snapshots show the adiabatic evolution of the σ(Cl−Cα) bond to the σ(Cα−F) bond. The bottom row shows the adiabatic evolution of the F−(2pz) orbital to the Cl−(3pz) orbital. The tiny orange and green dots represent the positions of F−and Cl−, respectively. Figure 3. Time evolution of the three chemically active natural orbitals for the syn-2 mechanism. The bottom line indicates the time stamps (fs) of the selected snapshots. Contour value = 0.1 au. The top row natural-orbitals’ snapshots show the adiabatic evolution of the σ(Cl−Cα) bond to the π(Cα−Cβ) bond. The middle row shows the adiabatic evolution of the F−(2pz) orbital to the Cl−(3pz) orbital, and the bottom row shows the adiabatic evolution of the σ(H−Cβ) bond to the σ(F−H) bond. The tiny orange and green dots represent the positions of F−and Cl−, respectively. See text for further details. The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter https://doi.org/10.1021/acs.jpclett.4c03010 J. Phys. Chem. Lett. 2024, 15, 12138−12143 12140 suggest that the syn-E2 mechanism resembles a two-step process where bond breaking and formation processes do not take place simultaneously but sequentially, opposite to SN2 (vide supra). anti-E2 Mechanism. The inspection of the chemically active natural orbitals of the anti-E2 mechanism, whose time evolution is shown in Figure 4, highlights not only some similarities with syn-E2 but also a number of profound and consequential differences, which will be discussed below. First, we observe that the reaction “starts” earlier than synE2. The interaction between the F−fluorine anion and the proton of σ(H−Cβ) is noticeable in the bottom row frame at t = 140 fs, whereas the onset for this interaction is 170 fs for synE2. Second, it is observed that the proton abstraction takes longer to complete. Namely, the “critical point” is arrived at t= 233 fs, which is to be compared with t= 194 fs for syn-E2. Further inspection of the evolution of the natural orbital describing the forming σ(H−F) bond (see bottom row frames from t= 130 fs to t= 233 fs) shows a rebound effect, meaning that such an σ(H−F) orbital reverses its development from t= 145 fs to t= 190 fs (see Figure S1 in the Supporting Information). Subsequently, from t= 190 fs onward the σ(H− F) bond formation resumes, and at t= 233 fs it gets fully formed and the Cβ−H is broken. This is very suggestive of the key role played by the H−Cβstretching vibrational modespecific excitation for this mechanism (vide infra). Third, similarly to the syn-E2 mechanism, the “critical point” marks the point for the evolution of the chemical bonding structure evolution to move on from the vicinity of Cβto Cα. From this critical point onward, the synchronous formation of the π(C−C) bond, at the expense of the delocalization of the Cβcarbanion’s electron pair over both carbon atoms, which necessarily turns out into a π-symmetry orbital due to the required orthogonality to the already existing σ-symmetry C− C bonding orbital, and the detachment of the Cl−chlorine anion, at the expense of the σ(Cl−Cα) bond breaking, take place in a short time of ∼6 fs. The difference is that now the “critical point” is arrived at t= 233 fs, ∼40 fs later than in synE2. This observation suggests that anti-E2 shows a late barrier (vide infra). Outline. Since the mechanism of a chemical reaction refers to the order in which chemical bonds are made and/or broken, 4 it follows naturally that the time evolution of the molecular natural orbitals is the proper “place” to look for chemically insightful descriptions of reaction mechanisms. The present study reports on the adiabatically relaxed molecular natural orbitals along selected trajectories for the canonical bimolecular (i) nucleophilic substitution and (ii) base-induced elimination reaction mechanisms of F−+ CH3CH2Cl, revealing how these analyses yield significant chemical insight into both competing reaction mechanisms. Our study takes advantage of the BO ab initio molecular dynamics where the adiabatic relaxation of the electronic structure is carried out by GNOF calculations. Thus, the nuclei are propagated following the classical equations of motion, and at each time step, the quantum mechanical electronic structure problem at the clamped nuclear positions is addressed using GNOF, which has been reported to be chemically accurate at describing molecular bond formation and breaking processes. Consequently, our procedure provides a detailed and reliable chemical description of the time evolution of the chemical bonding structure of the species involved in the reaction for each selected mechanism. We have shown that the SN2 nucleophilic substitution proceeds in one single step with the formation of the σ(F−Cα) chemical bond and the breaking of the σ(Cα−Cl) bond taking place synchronously. The chemical environment around the Cβ atom remains a spectator during the completion of the halogen substitution. The syn-E2 and the anti-E2 elimination mechanisms are substantially more complex because chemical action takes place on the two carbon atoms of CH3CH2Cl in a concerted but sequential way. The time-ordered chemical events include the abstraction of a proton from a σ(H−Cβ) bond by the incoming F−anion, and the subsequent formation of a transient carbanion at Cβ. During the completion of this first step, the Cαmoiety remains expectant. Next, the second step consists of the delocalization of the electron pair of the Cβ carbanion into the Cα, to form a π-type covalent bond between the two carbon atoms, and the Cl−detachment from the σ(Cl−Cα) bond. We observe, therefore, two sequential steps taking place in this order at the Cβand in the Cαatoms, which feature two decoupled bond breaking/formation processes occurring in different timeand real-space coordinate domains, neatly divided by the barrier epitomized by the hydrogen fluorine “stabilized” chlorinated carbanion transient structure. Furthermore, we observe that the proton abstraction by the incoming F−anion from the σ(H−Cβ) bond takes a longer Figure 4. Time evolution of the three chemically active natural orbitals for the anti-E2 mechanism. The bottom line indicates the time stamps (fs) of the selected snapshots. Contour value = 0.1 au. The top row natural-orbitals’ snapshots show the adiabatic evolution of the σ(Cl−Cα) bond to the π(Cα−Cβ) bond. The middle row shows the adiabatic evolution of the σ(H−Cβ) bond to the Cl−(3pz) orbital, and the bottom row shows the adiabatic evolution of the F−(2pz) orbital to the σ(F−H) bond. The tiny orange and green dots represent the positions of F−and Cl−, respectively. See text for further details. The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter https://doi.org/10.1021/acs.jpclett.4c03010 J. Phys. Chem. Lett. 2024, 15, 12138−12143 12141 time in anti-E2 (233 fs) than in syn-E2 (194 fs) and exhibits “rebound” effects in the former, in the sense that the proton moves back and forth between Cβand F−. This is very suggestive that vibrationally exciting the H−Cβstretching mode should enhance the proton abstraction and, eventually, should lead to anti-E2 dominance over both syn-E2 and SN2 mechanisms, as suggested elsewhere. 24,25 The second important implication of the longer time taken by the initial proton abstraction in anti-E2, is that the second step, i.e., the formation of the π(C−C) bond and detachment of Cl−, takes place at later times. Namely, the “barrier” for the transition from the first step to the second comes at later times in the anti-E2 bond relative to the SN2 mechanism. This is in accordance with the extension of the Polanyi rules 26 to polyatomic reactive processes. 25 Namely, increasing collision energy less efficiently activates late-barrier mechanisms than early barrier mechanisms and vice versa. Indeed, the experiments 15 of Meyer et al. convincingly show that for the herein studied reaction, anti-E2 ceases to be the most important mechanism at increased collision energies. ■APPENDIX Let us consider a mixed singlet state of N spin-paired electrons, specifically 44 in the F−+ CH3CH2Cl reaction. The GNOF is an electron-pairing based approach, 27 so the orbital space Ωis divided into N/2 mutually disjoint subspaces Ωg. Each orbital exclusively belongs to one Ωg, housing one strongly doubleoccupied orbital gand Ngweakly double-occupied orbitals. This partitioning ensures that the sum of the occupation numbers equals 2 for each subspace, with the trace of the oneparticle reduced density matrix equating the total number of electrons. The reconstruction functional for the two-particle reduced density matrix in terms of the occupation numbers leads to GNOF, as illustrated by the equation E E E E E el intra HF inter sta inter dyn inter = + + + (1) The intrapair component is formed by summing the energies of electron pairs, specifically, Ä Ç Å Å Å Å Å Å Å Å Å Å Å Å Å É Ö Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ Ñ E n H n n L2 ( , ) intra g N p p pp p q q p pq 1 /2 , g g = + = (2) where Hpp are the diagonal one-electron matrix elements of the kinetic energy and external potential operators, whereas Lpq = ⟨pp|qq⟩are the exchange-time-inversion integrals. 28 The matrix elements Π(nq,np) = c(nq)c(np), where c(np) is defined by the square root of the occupation numbers according to the following rule: l m o o o n o o o c n n p N n p N ( ) /2 /2 p p p = > (3) that is, the phase factor of c(np) is chosen to be +1 for the strongly occupied orbital of a given subspace Ωg, and −1 otherwise. The inter-subspace Hartree−Fock (HF) term is E n n J K(2 ) HF inter p q q p pq pq , = (4) where Jpq =⟨pq|pq⟩and Kpq =⟨pq|qp⟩are the Coulomb and exchange integrals, respectively. The prime in the summation indicates that only the inter-subspace terms are taken into account. The inter-subspace static component is written as E L(1 ) sta inter p q q p q p pq , b b = (5) where n h p p p = with the hole hp= 1−np. Finally, the intersubspace dynamic energy is E n n n n L(1 ) ( , ) dyn inter p q q p q d p d q d p d pq , b b = [ + ] (6) In eqs 5 and 6,Ωbdenotes the subspace composed of orbitals below the level N/2, whereas npdis the dynamic part of the occupation number npin accordance with the Pulay criterion that establishes an occupancy deviation of approximately 0.01 with respect to 1 or 0 for a natural orbital to contribute to the dynamic correlation. In the BO approximation, the total energy is cast as E=Eel + Enuc, where Eel represents the electronic energy calculated using (1), and Enuc is the nuclear energy E Z Z R( / ) A B A Bnuc AB =< . Here, ZAdenotes the atomic number of nucleus A, and RAB is the distance between nuclei Aand B. Considering that all natural orbitals are expanded in a fixed atomic basis set, the derivative of the total energy with respect to the coordinate xof nucleus Ais given by 29 E x E x H x D x S x d dA nuc A A A A = + + | (7) where Γμυ and Dμηυδ are the GNOF one-particle and twoparticle reduced density matrices, respectively, Sμυ =⟨μ|υ⟩is the overlap matrix, and λμυ are obtained from the reduced density matrices, all in the atomic orbital representation. The derivatives in eq 7 have an explicit dependence on the nuclear coordinate xA, so the force acting on each nucleus A(FA= −∇AE) can be obtained by a single static evaluation at each time step for the fixed nuclear positions at that instant. ■ASSOCIATED CONTENT * sı Supporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpclett.4c03010. SN2 reaction (MP4) syn-E2 reaction (MP4) anti-E2 reaction (MP4) Figure of interatomic distances versus time (PDF) Transparent Peer Review report available (PDF) ■AUTHOR INFORMATION Corresponding Author Mario Piris −Donostia International Physics Center (DIPC) &Kimika Fakultatea, Euskal Herriko Unibertsitatea (UPV/ EHU), 20018 Donostia, Euskadi, Spain; Basque Foundation for Science (IKERBASQUE), 48009 Bilbao, Euskadi, Spain; orcid.org/0000-0003-0222-2953; Email: mario.piris@ ehu.eus The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter https://doi.org/10.1021/acs.jpclett.4c03010 J. Phys. Chem. Lett. 2024, 15, 12138−12143 12142 Authors Xabier Lopez −Donostia International Physics Center (DIPC) &Kimika Fakultatea, Euskal Herriko Unibertsitatea (UPV/EHU), 20018 Donostia, Euskadi, Spain; orcid.org/0000-0002-2711-3588 Jesus M. Ugalde −Donostia International Physics Center (DIPC) &Kimika Fakultatea, Euskal Herriko Unibertsitatea (UPV/EHU), 20018 Donostia, Euskadi, Spain; orcid.org/0000-0001-8980-9751 Complete contact information is available at: https://pubs.acs.org/10.1021/acs.jpclett.4c03010 Notes The authors declare no competing financial interest. ■ACKNOWLEDGMENTS Financial support comes from the Eusko Jaurlaritza (Basque Government), ref.: IT1584-22 and from the Grant No. PID 2021-126714NB-I00, funded by MCIN/AEI/10.13039/ 501100011033. The authors are thankful for technical and human support provided by IZO-SGI SGIker of UPV/EHU and DIPC. ■REFERENCES (1) Steinfeld, J. I.; Fancisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics, 2nd ed.; Prentice Hall: Englewood Cliffs, NJ, 1999. (2) Tully, J. C. Molecular dynamics with electronic transitions. J. Chem. Phys. 1990,93, 1061−1071. (3) Paranjothy, M.; Sun, R.; Zhuang, Y.; Hase, W. L. Direct chemical dynamics simulations: coupling of classical and quasiclassical trajectories with electronic structure theory. WIREs Comput. Mol. Sci. 2013,3, 296−316. (4) Vallance, C. The fingerprints of reaction mechanisms. Nature 2017,546, 608−609. (5) Sherrill, C. Annu. Rep. Comput. Chem.; Elsevier, 2005; Vol. 1, Chapter 4, pp 45−56. (6) Piris, M. Global Natural Orbital Functional: Towards the Complete Description of the Electron Correlation. Phys. Rev. Lett. 2021,127, 233001. (7) Mitxelena, I.; Piris, M. Benchmarking GNOF against FCI in challenging systems in one, two, and three dimensions. J. Chem. Phys. 2022,156, 214102. (8) Mitxelena, I.; Piris, M. Assessing the Global Natural Orbital Functional Approximation on Model Systems with Strong Correlation. J. Chem. Phys. 2024,160, 204106. (9) Pan, H.; Liu, K.; Caracciolo, A.; Casavecchia, P. Crossed beam polyatomic reaction dynamics: Recent advances and new insights. Chem. Soc. Rev. 2017,46, 7517−7547. (10) Mugnai, M.; Cardini, G.; Schettino, V. Substitution and elimination reaction of F−with C2H5Cl: An ab initio molecular dynamics study. J. Phys. Chem. A 2003,107, 2540−2547. (11) Yang, L.; Zhang, J.; Xie, J.; Ma, X.; Zhang, L.; Zhao, C.; Hase, W. L. Competing E2 and SN2 Mechanisms for the F−+ CH3CH2I Reaction. J. Phys. Chem. A 2017,121, 1078−1085. (12) Carrascosa, E.; Meyer, J.; Zhang, J.; Stei, M.; Michaelsen, T.; Hase, W. L.; Yang, L.; Wester, R. Imaging dynamic fingerprints of competing E2 and SN2 reactions. Nat. Commun. 2017,8, 25. (13) Carrascosa, E.; Meyer, J.; Michaelsen, T.; Stei, M.; Wester, R. Conservation of direct dynamics in sterically hindered SN2/E2 reactions. Chem. Sci. 2018,9, 693−701. (14) Tasi, D. A.; Tokaji, C.; Czakó, G. A benchmark ab initio study of the complex potential energy surfaces of the OH−+ CH3CH2Y [Y = F, Cl, Br, I] reactions. Phys. Chem. Chem. Phys. 2021,23, 13526− 13534. (15) Meyer, J.; Tajti, V.; Carrascosa, E.; Gyori, T.; Stei, M.; Michaelsen, T.; Bastian, B.; Czako, G.; Wester, R. Atomistic dynamics of elimination and nucleophilic substitution disentangled for the F- + CH3CH2Cl reaction. Nat. Chem. 2021,13, 977−981. (16) Zhao, S.; Fu, G.; Zhen, W.; Wang, H.; Liu, M.; Yang, L.; Zhang, J. Nucleophile Effects on the E2/SN2 Competition for the X−+ CH3CH2Cl Reactions: A Theoretical Study. J. Phys. Chem. A 2023, 127, 3381−3389. (17) Zhao, S.; Wang, H.; Fu, G.; Zhen, W.; Liu, M.; Yang, L.; Zhang, J. E2/SN2 Selectivity Driven by Reaction Dynamics. Insight into Halogen Bonding. Precis. Chem. 2023,1, 507−515. (18) Nacsa, A. B.; Tokaji, C.; Czako, G. High-level analytical potential-energy-surface-based dynamics of the OH−+ CH3CH2Cl SN2 and E2 reactions in full (24) dimensions. Faraday Discuss. 2024, 251, 604−621. (19) Tajti, V.; Czakó, G. Benchmark ab Initio Characterization of the Complex Potential Energy Surface of the F−+ CH3CH2Cl Reaction. J. Phys. Chem. A 2017,121, 2847−2854. (20) Piris, M.; Mitxelena, I. DoNOF: an open-source implementation of natural-orbital-functional-based methods for quantum chemistry. Comput. Phys. Commun. 2021,259, 107651. Available at https://github.com/DoNOF/, documentation at https://donof. readthedocs.io/ (21) Beeman, D. Some multistep methods for use in molecular dynamics calculations. J. Comput. Phys. 1976,20, 130−139. (22) Dunning, T. H.; Dunning, T. H., Jr. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989,90, 1007−1023. (23) Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. III. The atoms aluminum through argon. J. Chem. Phys. 1993,98, 1358−1371. (24) Tajti, V.; Czakó, G. Vibrational mode-specific dynamics of the the F−+ CH3CH2Cl multi-channel reaction. Phys. Chem. Chem. Phys. 2022,24, 8166−8181. (25) Czakó, G.; Gruber, B.; Papp, D.; Tajti, V.; Tasi, D. A.; Yin, C. First-principles mode-specific reaction dynamics. Phys. Chem. Chem. Phys. 2024,26, 15818−15830. (26) Polanyi, J. C. Some Concepts in Reaction Dynamics. Science 1987,236, 680−690. (27) Piris, M. In Theoretical and Quantum Chemistry at the Dawn of the 21st Century; Chakraborty, T., Carbó-Dorca, R., Eds.; Apple Academic Press, 2018; pp 593−620. (28) Piris, M. A generalized self-consistent-field procedure in the improved BCS theory. J. Math. Chem. 1999,25, 47−54. (29) Mitxelena, I.; Piris, M. Analytic gradients for spin multiplets in natural orbital functional theory. J. Chem. Phys. 2020,153, 044101. ■NOTE ADDED AFTER ASAP PUBLICATION This paper was published ASAP on November 28, 2024, with an error in the caption to Figure 4. This was corrected in the version published ASAP on December 2, 2024. The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter https://doi.org/10.1021/acs.jpclett.4c03010 J. Phys. Chem. Lett. 2024, 15, 12138−12143 12143