Nonstatistical Reaction Dynamics.

Nonstatistical dynamics is important for many chemical reactions. The Rice-Ramsperger-Kassel-Marcus (RRKM) theory of unimolecular kinetics assumes a reactant molecule maintains a statistical microcanonical ensemble of vibrational states during its dissociation so that its unimolecular dynamics are time independent. Such dynamics results when the reactant's atomic motion is chaotic or irregular. Intrinsic non-RRKM dynamics occurs when part of the reactant's phase space consists of quasiperiodic/regular motion and a bottleneck exists, so that the unimolecular rate constant is time dependent. Nonrandom excitation of a molecule may result in short-time apparent non-RRKM dynamics. For rotational activation, the 2J + 1 K levels for a particular J may be highly mixed, making K an active degree of freedom, or K may be a good quantum number and an adiabatic degree of freedom. Nonstatistical dynamics is often important for bimolecular reactions and their intermediates and for product-energy partitioning of bimolecular and unimolecular reactions. Post-transition state dynamics is often highly complex and nonstatistical.

Mechanism and kinetics for the reaction of methyl peroxy radical with O2.

Quantum chemical calculations and dynamics simulations were performed to study the reaction between methyl peroxy radical (CH3O2) and O2. The reaction proceeds through three different pathways (1) H-atom abstraction, (2) O2 addition and (3) concerted H-atom shift and O2 addition reactions. The concerted H-atom shift and O2 addition pathway is the most favourable reaction both kinetically and thermodynamically. The major product channel formed from these reactions is H2CO + OH + O2. Trajectory calculations also confirm that H2CO + OH + O2 is the main product channel. An estimated rate constant expression for this reaction from master equation calculations is 4.20 × 1013 e-8676/T cm3 mole-1 s-1.

Direct Dynamics Simulations of the 3CH2 + 3O2 Reaction at High Temperature.

Direct dynamics simulations with the M06/6-311++G(d,p) level of theory were performed to study the 3CH2 + 3O2 reaction at 1000 K temperature on the ground state singlet surface. The reaction is complex with formation of many different product channels in highly exothermic reactions. CO, CO2, H2O, OH, H2, O, H, and HCO are the products formed from the reaction. The total simulation rate constant for the reaction at 1000 K is (1.2 ± 0.3) × 10-12 cm3 molecule-1 s-1, while the simulation rate constant at 300 K is (0.96 ± 0.28) × 10-12 cm3 molecule-1 s-1. The simulated product yields show that CO is the dominant product and the CO:CO2 ratio is 5.3:1, in good comparison with the experimental ratio of 4.3:1 at 1000 K. On comparing the product yields for the 300 and 1000 K simulations, we observed that, except for CO and H2O, the yields of the other products at 1000 K are lower at 300 K, showing a negative temperature dependence.

Sampling initial positions and momenta for nuclear trajectories from quantum mechanical distributions.

We compare algorithms to sample initial positions and momenta of a molecular system for classical trajectory simulations. We aim at reproducing the phase space quantum distribution for a vibrational eigenstate, as in Wigner theory. Moreover, we address the issue of controlling the total energy and the energy partition among the vibrational modes. In fact, Wigner's energy distributions are very broad, quite at variance with quantum eigenenergies. Many molecular processes depend sharply on the available energy, so a better energy definition is important. Two approaches are introduced and tested: the first consists in constraining the total energy of each trajectory to equal the quantum eigenenergy. The second approach modifies the phase space distribution so as to reduce the deviation of the single mode energies from the correct quantum values. A combination of the two approaches is also presented.

Collisional dynamics simulations revealing fragmentation properties of Zn(ii)-bound poly-peptide.

Chemical dynamics simulations are performed to study the collision induced gas phase unimolecular fragmentation of a model peptide with the sequence acetyl-His1-Cys2-Gly3-Pro4-Tyr5-His6-Cys7 (analogue methanobactin peptide-5, amb5) and in particular to explore the role of zinc binding in reactivity. Fragmentation pathways, their mechanisms, and collision energy transfer are discussed. The probability distributions of the pathways are compared with the results of the experimental IM-MS, MS/MS spectrum and previous thermal simulations. Collisional activation gives both statistical and non-statistical fragmentation pathways with non-statistical shattering mechanisms accounting for a relevant percentage of reactive trajectories, becoming dominant at higher energies. The tetra-coordination of zinc changes qualitative and quantitative fragmentation, in particular the shattering. The collision energy threshold for the shattering mechanism was found to be 118.9 kcal mol-1 which is substantially higher than the statistical Arrhenius activation barrier of 35.8 kcal mol-1 identified previously during thermal simulations. This difference can be attributed to the tetra-coordinated zinc complex that hinders the availability of the sidechains to undergo direct collision with the Ar projectile.

Dynamics of Pyrene-Dimer Association and Ensuing Pyrene-Dimer Dissociation.

To address the possible role of pyrene dimers in soot, chemical dynamics simulations are reported to provide atomistic details for the process of collisional association of pyrene dimers and ensuing decomposition of pyrene dimers. The simulations are performed at 600, 900, 1200, 1600, and 2000 K temperatures (T) with different collisional impact parameters (b; 0-18 Å) using the all-atom optimized potentials for liquid simulations intermolecular force field. Corresponding to each b, ensembles of 1000 trajectories are computed up to a maximum time of 110 ps at each T. Microcanonical association rate constants for the pyrene-dimerization processes decrease with an increase in T. The ensuing dissociation of the pyrene dimers is statistical and could be well represented by the Rice-Ramsperger-Kassel-Marcus theory of unimolecular dissociation. Fits of the dissociation rate constants versus the harmonic Rice-Ramsperger-Kassel equation revealed that partial energy randomization occurs among the inter- and intramolecular vibrational modes during the dissociation of pyrene dimers, whereas rotational and translational modes play a significant role. Based on the low probability of association and short lifetime at 1600 (∼13.3 ps) and 2000 (∼12.8 ps) K, it is concluded that pyrene dimers are unlikely to play any major role in soot nucleation processes.

Direct Dynamics Simulations of the Unimolecular Decomposition of the Randomly Excited 1CH2O2 Criegee Intermediate. Comparison with 3CH2 + 3O2 Reaction Dynamics.

The 3CH2 + 3O2 reaction has a quite complex ground state singlet potential energy surface (PES). There are multiple minima and transition states before forming the 10 possible reaction products. A previous direct chemical dynamics simulation at the UM06/6-311++G(d,p) level of theory ( J. Phys. Chem. A 2019, 123, 4360-4369) found that reaction on this PES is predominantly direct without trapping in the potential minima. The first minima 3CH2 + 3O2 encounters is that for the 1CH2O2 Criegee intermediate and statistical theory assumes the reactive system is trapped in this intermediate with a lifetime given by Rice-Ramsperger-Kassel-Marcus (RRKM) theory. In the work presented here, a direct dynamics simulation is performed with the above UM06 theory, with the trajectories initialized in the 1CH2O2 intermediate with a random distribution of vibrational energy as assumed by RRKM theory. There are substantial differences between the dynamics for 1CH2O2 dissociation and 3CH2 + 3O2 reaction. For the former there are four product channels, while for the latter there are seven in agreement with experiment. Product energy partitioning for the two simulations are in overall good agreement for the CO2 + H2 and CO + H2O product channels, but in significant disagreement for the HCO + OH product channel. Though 1CH2O2 is excited randomly in accord with RRKM theory, its dissociation probability is biexponential and not exponential as assumed by RRKM. In addition, the 1CH2O2 dissociation dynamics follow non-intrinsic reaction coordinate (non-IRC) pathways. An important finding is that the nonstatistical dynamics for the 3CH2 + 3O2 reaction give results in agreement with experiment.

Time-Dependent Perspective for the Intramolecular Couplings of the N-H Stretches of Protonated Tryptophan.

Quasi-classical direct dynamics simulations, performed with the B3LYP-D3/cc-pVDZ electronic structure theory, are reported for vibrational relaxation of the three NH stretches of the -NH3+ group of protonated tryptophan (TrpH+), excited to the n = 1 local mode states. The intramolecular vibrational energy relaxation (IVR) rates determined for these states, from the simulations, are in good agreement with the experiment. In accordance with the experiment, IVR for the free NH stretch is slowest, with faster IVR for the remaining two NH stretches which have intermolecular couplings with an O atom and a benzenoid ring. For the free NH and the NH coupled to the benzenoid ring, there are beats (i.e., recurrences) in their relaxations versus time. For the free NH stretch, 50% of the population remained in n = 1 when the trajectories were terminated at 0.4 ps. IVR for the free NH stretch is substantially slower than for the CH stretch in benzene. The agreement found in this study between quasi-classical direct dynamics simulations and experiments indicates the possible applicability of this simulation method to larger biological molecules. Because IVR can drive or inhibit reactions, calculations of IVR time scales are of interest, for example, in unimolecular reactions, mode-specific chemistry, and many photochemical processes.

Theoretical Study of the Dynamics of the HBr+ + CO2 → HOCO+ + Br Reaction.

The dynamics of the HBr+ + CO2 → HOCO+ + Br reaction was recently investigated with guided ion beam experiments under various excitations (collision energy of the reactants, rotational and spin-orbital states of HBr+, etc.), and their impacts were probed through the change of the cross section of the reaction. The potential energy profile of this reaction has also been accurately characterized by high-level ab initio methods such as CCSD(T)/CBS, and the UMP2/cc-pVDZ/lanl08d has been identified as an ideal method to study its dynamics. This manuscript reports the first ab initio molecular dynamics simulations of this reaction at two different collision energies, 8.1 kcal/mol and 19.6 kcal/mol. The cross sections measured from the simulations agree very well with the experiments measured with HBr+ in the 2∏1/2 state. The simulations reveal three distinct mechanisms at both collision energies: direct rebound (DR), direct stripping (DS), and indirect (Ind) mechanisms. DS and Ind make up 97% of the total reaction. The dynamics of this reaction is also compared with nucleophilic substitution (SN2) reactions of X- + CH3Y → CH3X + Y- type. In summary, this research has revealed interesting dynamics of the HBr+ + CO2 → HOCO+ + Br reaction at different collision energies and has laid a solid foundation for using this reaction to probe the impact of rotational excitation of ion-molecule reactions in general.

Exploring reactivity and product formation in N(4S) collisions with pristine and defected graphene with direct dynamics simulations.

Atomic nitrogen is formed in the high-temperature shock layer of hypersonic vehicles and contributes to the ablation of their thermal protection systems (TPSs). To gain atomic-level understanding of the ablation of carbon-based TPS, collisions of hyperthermal atomic nitrogen on representative carbon surfaces have recently be investigated using molecular beams. In this work, we report direct dynamics simulations of atomic-nitrogen [N(4S)] collisions with pristine, defected, and oxidized graphene. Apart from non-reactive scattering of nitrogen atoms, various forms of nitridation of graphene were observed in our simulations. Furthermore, a number of gaseous molecules, including the experimentally observed CN molecule, have been found to desorb as a result of N-atom bombardment. These results provide a foundation for understanding the molecular beam experiment and for modeling the ablation of carbon-based TPSs and for future improvement of their properties.

Comparison of intermolecular energy transfer from vibrationally excited benzene in mixed nitrogen-benzene baths at 140 K and 300 K.

Gas phase intermolecular energy transfer (IET) is a fundamental component of accurately explaining the behavior of gas phase systems in which the internal energy of particular modes of molecules is greatly out of equilibrium. In this work, chemical dynamics simulations of mixed benzene/N2 baths with one highly vibrationally excited benzene molecule (Bz*) are compared to experimental results at 140 K. Two mixed bath models are considered. In one, the bath consists of 190 N2 and 10 Bz, whereas in the other bath, 396 N2 and 4 Bz are utilized. The results are compared to results from 300 K simulations and experiments, revealing that Bz*-Bz vibration-vibration IET efficiency increased at low temperatures consistent with longer lived "chattering" collisions at lower temperatures. In the simulations, at the Bz* excitation energy of 150 kcal/mol, the averaged energy transferred per collision, ⟨ΔEc⟩, for Bz*-Bz collisions is found to be ∼2.4 times larger in 140 K than in 300 K bath, whereas this value is ∼1.3 times lower for Bz*-N2 collisions. The overall ⟨ΔEc⟩, for all collisions, is found to be almost two times larger at 140 K compared to the one obtained from the 300 K bath. Such an enhancement of IET efficiency at 140 K is qualitatively consistent with the experimental observation. However, the possible reasons for not attaining a quantitative agreement are discussed. These results imply that the bath temperature and molecular composition as well as the magnitude of vibrational energy of a highly vibrationally excited molecule can shift the overall timescale of rethermalization.

The generality of the GUGA MRCI approach in COLUMBUS for treating complex quantum chemistry.

The core part of the program system COLUMBUS allows highly efficient calculations using variational multireference (MR) methods in the framework of configuration interaction with single and double excitations (MR-CISD) and averaged quadratic coupled-cluster calculations (MR-AQCC), based on uncontracted sets of configurations and the graphical unitary group approach (GUGA). The availability of analytic MR-CISD and MR-AQCC energy gradients and analytic nonadiabatic couplings for MR-CISD enables exciting applications including, e.g., investigations of π-conjugated biradicaloid compounds, calculations of multitudes of excited states, development of diabatization procedures, and furnishing the electronic structure information for on-the-fly surface nonadiabatic dynamics. With fully variational uncontracted spin-orbit MRCI, COLUMBUS provides a unique possibility of performing high-level calculations on compounds containing heavy atoms up to lanthanides and actinides. Crucial for carrying out all of these calculations effectively is the availability of an efficient parallel code for the CI step. Configuration spaces of several billion in size now can be treated quite routinely on standard parallel computer clusters. Emerging developments in COLUMBUS, including the all configuration mean energy multiconfiguration self-consistent field method and the graphically contracted function method, promise to allow practically unlimited configuration space dimensions. Spin density based on the GUGA approach, analytic spin-orbit energy gradients, possibilities for local electron correlation MR calculations, development of general interfaces for nonadiabatic dynamics, and MRCI linear vibronic coupling models conclude this overview.

Addressing an instability in unrestricted density functional theory direct dynamics simulations.

In Density Functional Theory (DFT) direct dynamics simulations with Unrestricted Hartree Fock (UHF) theory, triplet instability often emerges when numerically integrating a classical trajectory. A broken symmetry initial guess for the wave function is often used to obtain the unrestricted DFT potential energy surface (PES), but this is found to be often insufficient for direct dynamics simulations. An algorithm is described for obtaining smooth transitions between the open-shell and the closed-shell regions of the unrestricted PES, and thus stable trajectories, for direct dynamics simulations of dioxetane and its •OCH2 -CH2 O• singlet diradical. © 2018 Wiley Periodicals, Inc.

Pronounced changes in atomistic mechanisms for the Cl- + CH3I SN2 reaction with increasing collision energy.

In a previous direct dynamics simulation of the Cl- + CH3I → ClCH3 + I- SN2 reaction, predominantly indirect and direct reaction was found at collision energies Erel of 0.20 and 0.39 eV, respectively. For the work presented here, these simulations were extended by studying the reaction dynamics from Erel of 0.15 to 0.40 eV in 0.05 eV intervals. A transition from a predominantly indirect to direct reaction is found for Erel of 0.27-0.28 eV, a finding consistent with experiment. The simulation results corroborate the understanding that in experiments indirect reaction is characterized by small product translational energies and isotropic scattering, while direct reaction has higher translational energies and anisotropic scattering. The traditional statistical theoretical model for the Cl- + CH3I SN2 reaction assumes the Cl--CH3I pre-reaction complex (A) is formed, followed by barrier crossing, and then formation of the ClCH3-I- post-reaction complex (B). This mechanism is seen in the dynamics, but the complete atomistic dynamics are much more complex. Atomistic SN2 mechanisms contain A and B, but other dynamical events consisting of barrier recrossing (br) and the roundabout (Ra), in which the CH3-moiety rotates around the heavy I-atom, are also observed. The two most important mechanisms are only formation of A and Ra + A. The simulation results are compared with simulations and experiments for Cl- + CH3Cl, Cl- + CH3Br, F- + CH3I, and OH- + CH3I.

Unimolecular Rate Constants versus Energy and Pressure as a Convolution of Unimolecular Lifetime and Collisional Deactivation Probabilities. Analyses of Intrinsic Non-RRKM Dynamics.

Following work by Slater and Bunker, the unimolecular rate constant versus collision frequency, kuni(ω, E), is expressed as a convolution of unimolecular lifetime and collisional deactivation probabilities. This allows incorporation of nonexponential, intrinsically non-RRKM, populations of dissociating molecules versus time, N( t)/ N(0), in the expression for kuni(ω, E). Previous work using this approach is reviewed. In the work presented here, the biexponential f1 exp(- k1 t) + f2 exp(- k2 t) is used to represent N( t)/ N(0), where f1 k1 + f2 k2 equals the RRKM rate constant k( E) and f1 + f2 = 1. With these two constraints, there are two adjustable parameters in the biexponential N( t)/ N(0) to represent intrinsic non-RRKM dynamics. The rate constant k1 is larger than k( E) and k2 is smaller. This biexponential gives kuni(ω, E) rate constants that are lower than the RRKM prediction, except at the high and low pressure limits. The deviation from the RRKM prediction increases as f1 is made smaller and k1 made larger. Of considerable interest is the finding that, if the collision frequency ω for the RRKM plot of kuni(ω, E) versus ω is multiplied by an energy transfer efficiency factor ßc, the RRKM kuni(ω, E) versus ω plot may be scaled to match those for the intrinsic non-RRKM, biexponential N( t)/ N(0), plots. This analysis identifies the importance of determining accurate collisional intermolecular energy transfer (IET) efficiencies.

Direct Dynamics Simulations of the CH2 + O2 Reaction on the Ground- and Excited-State Singlet Surfaces.

In a previous work [ Lakshmanan , S. ; J. Phys. Chem. A 2018 , 122 , 4808 - 4818 ], direct dynamics simulations at the M06/6-311++G(d,p) level of theory were reported for 3CH2 (X3B1) + 3O2 (X3∑g-) reaction on its ground-state singlet potential energy surface (PES) at 300 K. However, further analyses revealed the simulations are unstable for the 3CH2 (X3B1) + 3O2 (X3∑g-) reactants on the ground-state singlet surface and the trajectories reverted to an excited-state singlet surface for the 1CH2 (ã1A1) + 1O2 (b1∑g+) reactants. Thus, the dynamics reported previously are for this excited-state singlet PES. The PESs for the 3CH2 (X3B1) + 3O2 (X3∑g-) and 1CH2 (ã1A1) + 1O2 (b1∑g+) reactants are quite similar, and this provided a means to perform simulations for the 3CH2 (X3B1) + 3O2 (X3∑g-) reactants on the ground-state singlet PES at 300 K, which are reported here. The reaction dynamics are quite complex with seven different reaction pathways and nine different products. A consistent set of product yields have not been determined experimentally, but the simulation yields for the H atom, CO, and CO2 are somewhat lower, higher, and lower respectively, than the recommended values. The yields for the remaining six products agree with experimental values. Product decomposition was included in determining the product yields. The simulation 3CH2 + 3O2 rate constant at 300 K is only 3.4 times smaller than the recommended value, which may be accommodated if the 3CH2 + 3O2 → 1CH2O2 potential energy curve is only 0.75 kcal/mol more attractive at the variational transition state for 3CH2 + 3O2 → 1CH2O2 association. The simulation kinetics and dynamics for the 3CH2 + 3O2 and 1CH2 + 1O2 reactions are quite similar. Their rate constants are statistically the same, an expected result, since their transition states leading to products have energies lower than that of the reactants and the attractive potential energy curves for 3CH2 + 3O2 → 1CH2O2 and 1CH2 + 1O2 → 1CH2O2 are nearly identical. The product yields for the 3CH2 + 3O2 and 1CH2 + 1O2 reactions are also nearly identical, only differing for the CO2 yield. The reaction dynamics on both surfaces are predominantly direct, with negligible trapping in potential energy minima, which may be an important contributor to their nearly identical product yields.

Chemical Dynamics Simulation of Energy Transfer: Propylbenzene Cation and N2 Collisions.

Collisional energy transfer of highly vibrationally excited propylbenzene cation in a N2 bath has been studied with chemical dynamics simulations. In this work, an intermolecular potential of propylbenzene cation interacting with N2 was developed from SCS-MP2/6-311++G** ab initio calculations. Using a particle swarm optimization algorithm, the ab initio results were simultaneously fit to a sum of three two-body potentials, consisting of C a-N, C b-N, and H-N, where C a is carbon on the benzene ring and C b is carbon on the propyl side chain. Using the developed intermolecular potential, classical trajectory calculations were performed with a 100.1 kcal/mol excitation energy at 473 K to compare with experiment. Varying the density of the N2 bath, the single collision limit of propylbenzene cation with the N2 bath was obtained at a density of 20 kg/m3 (28 atm). For the experimental excitation energy and in the single collision limit, the average energy transferred per collision, ⟨Δ Ec⟩, is 1.04 ± 0.04 kcal/mol and in good agreement with the experimental value of 0.82 kcal/mol.

Direct Dynamics Simulations of Fragmentation of a Zn(II)-2Cys-2His Oligopeptide. Comparison with Mass Spectrometry Collision-Induced Dissociation.

Abnormalities in zinc metabolism have been linked to many diseases, including different kinds of cancers and neurological diseases. The present study investigates the fragmentation pathways of a zinc chaperon using a model peptide with the sequence acetyl-His1-Cys2-Gly3-Pro4-Tyr5-His6-Cys7 (analog methanobactin peptide-5, amb5). DFT/M05-2X and B3LYP geometry optimizations of [amb5-3H+Zn(II)]- predicted three lowest energy conformers with different chelating motifs. Direct dynamics simulations, using the PM7 semiempirical electronic structure method, were performed for these conformers, labeled a, b, and c, to obtain their fragmentation pathways at different temperatures in the range 1600-2250 K. The simulation results were compared with negative ion mode mass spectrometry experiments. For conformer a, the number of primary dissociation pathways are 11, 14, 24, 70, and 71 at 1600, 1750, 1875, 2000, and 2250 K, respectively. However, there are only 6, 10, 13, 14, and 19 pathways corresponding to these temperatures that have a probability of 2% or more. For conformer b, there are 67 pathways at 2000 K and 71 pathways at 2250 K. For conformer c, 17 pathways were observed at 2000 K. For conformer a, for two of the most common pathways involving C-C bond dissociation, Arrhenius parameters were calculated. The frequency factors and activation energies are smaller than those for C-C homolytic dissociation in alkanes due to increased stability of the product ions as a result of hydrogen bonding. The activation energies agree with the PM7 barriers for the C-C dissociations. Comparison of the simulation and experimental fragmentation ion yields shows the simulations predict double or triple cleavages of the backbone with Zn(II) retaining its binding sites, whereas the experiment exhibits single cleavages of the backbone accompanied by cleavage of two of the Zn(II) binding sites, resulting in b- and y-type ions.

Potential Energy Curves for Formation of the CH2O2 Criegee Intermediate on the 3CH2 + 3O2 Singlet and Triplet Potential Energy Surfaces.

The potential energy curves (PECs) for the interaction of 3CH2 with 3O2 in singlet and triplet potential energy surfaces (PESs) leading to singlet and triplet Criegee intermediates (CH2OO) are studied using electronic structure calculations. The bonding mechanism is interpreted by analyzing the ground state multireference configuration interaction (MRCI) wave function of the reacting species and at all points along the PES. The interaction of 3CH2 with 3O2 on the singlet surface leads to a flat long-range attractive PEC lacking any maxima or minima along the curve. The triplet surface stems into a maximum along the curve resulting in a transition state with an energy barrier of 5.3 kcal/mol at CASSCF(4,4)/cc-pVTZ level. The resulting 3CH2OO is less stable than the 1CH2OO. In this study, the biradical character (ß) is used as a measure to understand the difference in the topology of the singlet and triplet PECs and the relation of the biradical nature of the species with their structures. The 3CH2OO has a larger biradical character than 1CH2OO, and because of the larger bond order of 1CH2OO, the C-O covalent bond becomes harder to break, thereby stabilizing 1CH2OO. Thus, this study provides insights into the shape of the PEC obtained from the reaction between 3CH2 and 3O2 in terms of their bonding nature and from the shape of the curves, the temperature dependence or independence of the rate of the reaction is discussed.

l-Cysteine Modified by S-Sulfation: Consequence on Fragmentation Processes Elucidated by Tandem Mass Spectrometry and Chemical Dynamics Simulations.

Low-energy collision-induced dissociation (CID) of deprotonated l-cysteine S-sulfate, [cysS-SO3]-, delivered in the gas phase by electrospray ionization, has been found to provide a means to form deprotonated l-cysteine sulfenic acid, which is a fleeting intermediate in biological media. The reaction mechanism underlying this process is the focus of the present contribution. At the same time, other novel species are formed, which were not observed in previous experiments. To understand fragmentation pathways of [cysS-SO3]-, reactive chemical dynamics simulations coupled with a novel algorithm for automatic determination of intermediates and transition states were performed. This approach has allowed the identification of the mechanisms involved and explained the experimental fragmentation pathways. Chemical dynamics simulations have shown that a roaming-like mechanism can be at the origin of l-cysteine sulfenic acid.