Next Article in Journal
The Prebiotic Activity of Simulated Gastric and Intestinal Digesta of Polysaccharides from the Hericium erinaceus
Next Article in Special Issue
First-Principles Study of the Reaction between Fluorinated Graphene and Ethylenediamine
Previous Article in Journal
GC-MS Characterization of Volatile Flavor Compounds in Stinky Tofu Brine by Optimization of Headspace Solid-Phase Microextraction Conditions
Previous Article in Special Issue
Reaction between Indazole and Pd-Bound Isocyanides—A Theoretical Mechanistic Study
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

A Trajectory-Based Method to Explore Reaction Mechanisms

by
Saulo A. Vázquez
1,
Xose L. Otero
2 and
Emilio Martinez-Nunez
1,*
1
Departamento de Química Física, Facultade de Química, Campus Vida, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain
2
Unidade de Bioestadística, Facultade de Medicina, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Spain
*
Author to whom correspondence should be addressed.
Molecules 2018, 23(12), 3156; https://doi.org/10.3390/molecules23123156
Submission received: 27 October 2018 / Revised: 23 November 2018 / Accepted: 29 November 2018 / Published: 30 November 2018
(This article belongs to the Special Issue Theoretical Investigations of Reaction Mechanisms)

Abstract

:
The tsscds method, recently developed in our group, discovers chemical reaction mechanisms with minimal human intervention. It employs accelerated molecular dynamics, spectral graph theory, statistical rate theory and stochastic simulations to uncover chemical reaction paths and to solve the kinetics at the experimental conditions. In the present review, its application to solve mechanistic/kinetics problems in different research areas will be presented. Examples will be given of reactions involved in photodissociation dynamics, mass spectrometry, combustion chemistry and organometallic catalysis. Some planned improvements will also be described.

Graphical Abstract

1. Introduction

Theoretical studies of reaction mechanisms can greatly benefit nowadays by leveraging the surge of automated methods developed in the last few years [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58]. The idea of these new computational protocols is to substitute human intervention by less error-prone and less tedious automated algorithms. The automated methodologies range from chemical heuristics to the use of artificial forces to boost chemical reactions.
Our group has contributed with the development of a method called tsscds [43,44,45,46,47], which is based on accelerated molecular dynamics (MD), as are some others [29,30]. In our trajectories, the bonds of the molecule(s) are broken/formed thanks to large amounts of energy placed in each normal mode/atom of the system [45]. The distinctive feature of tsscds compared to others is the primary target of the post-processing analysis: the search for transition states (TS) rather than minima. Additionally, having determined the TS of a given process, its rate can easily be determined using transition state theory (TST) [59,60,61,62]. Thus, finding the relevant TSs on a given potential energy surface (PES), as our method does, is a subject of fundamental importance in chemistry.
In tsscds, after completion of a trajectory, an algorithm named bond breaking/formation search (BBFS) [45] is employed to select good TS guess structures, which are then optimized using Eigenvector Following (EF) [63]. In particular, the adjacency matrix, which indicates whether pairs of atoms form a bond, is monitored along each trajectory to identify the atoms/bonds involved in all chemical reactions taking place. Then, for each of the selected candidates, a partial optimization is firstly carried out by freezing the atoms involved in the reaction. The partially-optimized structure is subsequently subjected to TS optimization using the EF algorithm. The resulting TSs are then connected with the minima using intrinsic reaction coordinate (IRC) calculations [64]. Finally, tsscds also features a Kinetic Monte Carlo [65] module that provides the desired kinetic information using the network of TSs and minima. The source code can be downloaded from: http://forge.cesga.es/wiki/g/tsscds/HomePage.
The method has been successfully employed to study reactions involved in combustion [66,67], photolysis [68,69,70], mass spectrometry [71] and organometallic catalysis [43]. The aim of this review is to go over several examples where tsscds is employed to either discover new mechanisms and/or to explain the experiments. For detailed comparisons among different methods for exploring reaction space, the reader is referred to two recent reviews [58,72]. Additionally, in the last section, some planned improvements to enhance the efficiency/efficacy or to expand the scope of tsscds will be described.

2. Method

The method tsscds has been recently put forward by one of the authors as an automated tool to discover reaction mechanisms [44,45]. The basic idea behind tsscds is to run accelerated MD simulations with the aim to break/form bonds within a few hundred femtoseconds. The simulations are called “accelerated” because the molecules experience breakage or formation of new bonds very rapidly thanks to large amounts of vibrational energy placed in each normal mode of the system. In particular, a range of vibrational energies of ~20–50 kcal/mol per normal mode is initially employed. However, this range is automatically adjusted to attain at least 60% reactive trajectories in the MD simulations. Although the default option is to excite all vibrational modes of the system (using microcanonical normal mode sampling [73]), the user can decide to heat only one part of the system selecting a few normal modes to be initially excited. The latter option can be particularly useful for large systems.
The trajectory results are then analyzed with a post-processing algorithm (named BBFS), which identifies geometries with partly formed/broken bonds. Those structures serve as TS candidates in subsequent transition state optimizations. As detailed below, BBFS is based on the adjacency matrix, a Graph Theory object that has been employed in other successful automated methods like the one developed by Zimmerman [16]. Similar ideas have also been recently employed to analyze changes in conformations occurring in MD simulations [74].
Once the TSs are optimized, a reaction network is constructed by computing the intrinsic reaction coordinates (IRCs) [64] connecting TSs with intermediates [64]. The method employs two levels of theory: semi-empirical and ab initio/DFT. The semi-empirical calculations are performed to run the MD simulations and to obtain approximate TSs structures, while a higher level of theory is used to re-optimize the TSs and run IRC calculations. Two different electronic structure programs are employed: MOPAC2016 [75] and Gaussian09 [76] for the semi-empirical and ab initio/DFT calculations, respectively.
Unlike other automated methods like GRRM [42], our methodology has been employed so far to study only the ground electronic state. This is in part due to the fact that, currently, the potential energy and gradients can only be calculated at the semiempirical level of theory. The following is a description of the graph-theoretic tools and kinetic models employed in our method.

2.1. Graph Theory

A number of graph theoretic tools are employed at various stages of the procedure to find transition states (TS), screen their structures and construct a reaction network. Specifically, the time dependence of the adjacency matrix A is employed to discriminate TS-like geometries along the trajectories. The elements of this matrix are defined as:
a i j = { 1   if   r i j < r i j ref 0     otherwise
with r i j being the distance between atoms i and j , and r i j ref a reference value that sets the upper limit for the bond length between the pair; in practice r i j ref is taken 20% greater than the sum of the covalent radii of i and j [45]. Thus, for an N - atom system, A is a N × N symmetric matrix with zeros on its diagonal.
Additionally, a weighted adjacency matrix A w is also employed in tsscds, whose off-diagonal elements are defined as:
a i j w = 1 ( r i j / r i j ref ) n 1 ( r i j / r i j ref ) m
Values of 6 and 12 have been employed in previous work for n and m , respectively [44]. Matrix A w contains information on the 3D geometry of the molecule [77] and its eigenvalues and eigenvectors can be employed to construct the so-called SPRINT coordinates [77]. An important property of these coordinates is their invariance with respect to translation, rotation and permutation of atoms, which makes them good molecular descriptors in trajectory-based methods. SPRINT coordinates are employed in tsscds to remove redundant structures.
Another matrix employed to determine the number of fragments in the system is the Laplacian, which is defined as:
L ( w ) = D A ( w )
where D is the so-called degree matrix [44], whose elements are defined as:
d i j = { deg ( v i ) if   i = j 0     otherwise    
where the degree deg ( v i ) of an atom counts the number of contacts. The superscript ( w ) on L and A indicates that the corresponding matrix can either be weighted or not. For a non-weighted graph, the lowest eigenvalue of the Laplacian λ 1 is always zero, and the total number of zero eigenvalues determines the number of fragments of the system. For a weighted graph, an upper threshold for λ 1 w is employed to identify fragmented structures [44]. The smallest non-zero eigenvalue is called the spectral gap, which is a measure of the degree of fragmentation of the structure. Thus, a small value of the spectral gap is associated with structures presenting non-covalent bonds (like van der Waals complexes), which are usually of no interest in chemical dynamics and kinetics.
The invariance of the SPRINT coordinates upon atom permutation is very important for the analyses of trajectories, where scrambling of atoms is frequent, as stated above. However, since the identity of each atom is absent in the adjacency matrix, SPRINT coordinates are identical for two structures where two non-equivalent atoms swap positions. For that reason, another type of molecular descriptor, based on a modified (weighted or not) adjacency matrix, is employed in tsscds. This new matrix, denoted as A Z ( w ) , contains the atomic numbers Z i of the atoms on the diagonal:
a Z , i j ( w ) = { a i j ( w )   if   i j 1 + Z i 10     if   i = j
The expression for the diagonal elements is chosen to provide values comparable to the off-diagonal ones. Most importantly, the eigenvalues of this new matrix are only invariant with respect to the permutation of like atoms, and it is widely employed in tsscds.

2.2. Kinetics Simulations

The kinetics module of tsscds calculates rate constants for all the elementary steps and solves the set of first-order differential equations that describe the time evolution of all species (usually known as chemical master equation).
The rate constants can either be obtained as a function of temperature or energy. In the former case, transition state theory is employed [59,60,61,62]:
k ( T ) = σ k B T h ( R T p 0 ) Δ n e Δ G R T
where σ is the reaction path degeneracy, T is the temperature, h is Planck’s constant, Δ G is the free energy of activation, p 0 is 1 bar and Δ n = 1 (0) for bimolecular (unimolecular) reactions. The reaction path degeneracy is calculated as σ = m T S m , where m and m T S are the number of optical isomers of the reactant and transition states, respectively [78].
By contrast, the microcanonical rate constants are computed according to RRKM theory [78]:
k ( E ) = σ W T S ( E ) h ρ ( E )
where W T S ( E ) is the sum of states at the TS, ρ ( E ) is the density of states at the reactant, and E is the excitation energy of the system. The sums and densities of states are evaluated by direct count of the harmonic vibrational states using the Beyer-Swinehart algorithm. Once all state-to-state rates are determined, the chemical master equation is solved using Kinetic Monte Carlo simulations [65].

3. Overview of the Applications of Tsscds

The tsscds methodology has been employed in our lab to elucidate reaction mechanisms involved in photodissociation dynamics, mass spectrometry, combustion and organometallic catalysis, and in this section, several examples of each type are reviewed.

3.1. Photodissociation Dynamics

The dissociation of molecules can be promoted by using a laser source, which is known as photodissociation. Although many photodissociations take place in excited states, important mechanisms may occur in the ground electronic state following internal conversion. One of the quantities of interest is the product yield, which is usually determined in the experiments. The understanding of the dissociation channels in organic compounds has greatly benefited from the interplay between photolysis experiments and computational studies [70,79,80,81,82,83,84,85,86,87,88,89,90,91,92].
In this section, we summarize the results obtained with our automated method for systems that have also been studied in photodissociation experiments, highlighting the most important conclusions. In particular, the dissociation channels of formaldehyde, formic acid, vinyl cyanide, acrolein, acryloyl chloride and methyl cyanoformate were studied with our tsscds methodology.
Formaldehyde was employed as a benchmark system to test tsscds. The system had been previously studied with other automated methods like the scaled hypersphere search [33] and the global reaction route mapping (GRRM) [35]. The results obtained with all algorithms are comparable, and the kinetically-relevant stationary points are found using any procedure.
The study of the dissociation channels of formic acid (CO2H2) with tsscds revealed the existence of a new TS for the water-gas shift reaction (WGSR: CO + H2O → CO2 + H2) [45]. By contrast, GRRM predicted three consecutive steps for the shortest path of the WGSR [35]. The discovery of the new TS is a consequence of the highly non-IRC [93] nature of the trajectories employed in tsscds [45]; in other words, IRC jumps are not uncommon events [94]. This exemplifies one of the advantages of using trajectory-based methods to discover new reactions: we are not restricted to unimolecular reactions and the only constrain to discover new processes is the molecular formula of the system. Additionally, the large amounts of vibrational energy put in the normal modes enhances configurational space sampling in tsscds, which permits the exploration of all types of reactions.
Our automated computational study on the dissociation of vinyl cyanide (VCN) [70] provides a HCN/HNC branching ratio in nearly perfect agreement with experiments for an excitation energy of 148 kcal/mol [95]. Besides the traditional 3-center and 4-center elimination mechanisms found in many HX eliminations from CH2=CHX systems, a new HCN elimination pathway involving three TSs was discovered in the tsscds study. The new mechanism involves three TSs and two intermediates and is shown in Figure 1.
Although alternative routes for HX elimination were also found for other ethylene analogues, those pathways involved high-energy TSs and were not competitive with the conventional 3-center and 4-center channels. This was the first time a new HX elimination channel competes with the well-known 3-center and 4-center processes in the dissociation of CH2=CHX species.
Figure 2 shows the product yields as a function of excitation energy obtained in our kinetic simulations from VCN. As seen in the figure, at low excitation energies (<150 kcal/mol) the new channel (red) is more important than the 4-center channel (green) and accounts for half of the HCN eliminations when the excitation energy is 110 kcal/mol.
The tsscds methodology was also employed to study the dissociation of acrolein (ACRL, C3H4O), which comprises many different fragmentation channels involving more than 250 transition states and 66 minima [44]. This system was studied with an enhanced procedure (now fully integrated in the method) consisting in the initialization of the MD simulations from multiple minima. In this new procedure the method works in an iterative manner. In the first iteration all MD simulations start from a starting structure, but once some TSs and intermediates are located, subsequent iterations utilize not only the starting equilibrium structure but also the newly generated intermediates to initialize the MD simulations. Compare to a single-minimum initialization, the use of multiple minima to start the dynamics ensures a better sampling of the PES of the system.
The potential energy surface of the C3H4O system is very complex and the 32 equilibrium structures (not including conformers) shown in Figure 3 were found with tsscds, with ACRL being the global minimum. To exemplify the importance of automated reaction discovery methods, we compare our results with those obtained by Chin et al. [96], who manually located equilibrium structures and TSs. Using the same levels of theory as in our study, Chin et al. only found 6 of the 66 minima obtained with tsscds. Most importantly, the relative product abundances obtained with tsscds at 148 kcal/mol (the energy corresponding to the experimental wavelength of 193 nm) are much closer to the experimental results than the computational results of Chin et al., as seen in Table 1.
Another system that attracted our attention was acryloyl chloride (AC). Overall, around 700 stationary points were found using our tsscds strategy. Of all possible dissociation channels from AC, experiments focus on the HCl dissociations. The use of our automated procedure led to the discovery of the three new HCl dissociation TSs [69] displayed in Figure 4; the figure also shows the AC equilibrium structure. The highest-energy TSs (TS2 and TS3) correspond to three-body dissociations leading to acetylene, carbon monoxide and hydrogen chloride, and they only become important at high excitation energies. By contrast, HCl elimination over TS1 is predominant at the experimental conditions (148 kcal/mol) [98], showing again that tsscds is capable of finding competitive pathways.
Finally, with the aim of exploring possible sources of HCN and HNC in astrophysical environments, the dissociation channels of methyl cyanoformate (MCF) were probed with tsscds, excited state calculations and photolysis experiments [68]. In particular, time-resolved infrared spectroscopy measurements indicate that both HCN and HNC are formed after the 193-nm photolysis of MCF [68]. The excited state calculations suggest that most of the dissociations take place in the S2 excited state leading to CH3O + NCCO via a Norrish type I reaction, in agreement with experiment. However, our calculations are also consistent with cascading internal conversion from S2 to produce vibrationally excited ground state MCF.
To study the dissociation of vibrationally excited MCF molecules in the S0 electronic state, tsscds was employed. Our approach assumes that, after the internal conversion process, intramolecular vibrational redistribution is fast enough to ensure RRKM behavior. With the tsscds procedure several HNC and HCN mechanisms are found, and Figure 5 shows the kinetically-relevant ones at 148 kcal/mol. The kinetic simulations predict a HNC/HCN branching ratio of 0.01, which is in semiquantitative agreement with that determined in the experiments (≈0.07). The work provides further insights into the intriguing observation of overabundance of HNC in astrophysical environments.

3.2. Mass Spectrometry

The prediction of mass spectra remains much of a challenge for the community of computational chemists. The common computational approaches employed for such endeavor include statistical rate theory calculations, MD simulations and electronic structure calculations [99,100,101,102,103,104,105,106,107,108,109,110,111,112,113]. Our automated method is very useful in this regard and can easily be coupled with MD simulations of collisions to generate theoretically-based mass spectra as described below.
In particular, tsscds was employed to simulate mass spectrometry (MS) experiments of protonated uracil, [uracil]H+. Our computational results indicate that the decomposition of [uracil]H+ involves more than one thousand stationary points and 751 elementary reactions [71]. Branching ratios for the different fragmentation channels can be automatically obtained from tsscds. However, these fractions are a function of the ion’s internal energy and cannot be directly compared with MS experiments, where the collision energy in the center-of-mass framework ( E c o m ) is employed instead. For that reason the tsscds results were combined with collisional dynamics simulations [71], which provide the fraction of E c o m transferred to the ion’s internal energy.
The resulting computationally-predicted product abundances (dashed lines) are compared in Figure 6 with the experimental ones (solid lines). As seen in the figure, for the predominant dissociation channels, the computationally-predicted product abundances are in qualitative agreement with experiment, and formation of HNCO (black), NH3 (red), H2O (green) and HNCOH+ (blue) are the major channels. Discrepancies with experiment can be attributed to the possible existence of well-known non-statistical behavior in many collision-induced dissociations [100,114], which cannot be captured with our statistical model.

3.3. Combustion Chemistry

Modeling the combustion reactions of oxygenated fuels is of great interest due to their potential use as alternatives to conventional petroleum-based fuels. To investigate combustion mechanisms, it is important to use kinetic models and perform computer simulations as a complement to experimental determinations, due to the tremendous complexity of these chemical processes. In general, different approximations are employed in combustion simulations to handle the complicated mechanisms. One of these simplifications consist of considering only the lowest energy rotamers of the involved species, which can lead to large errors in the calculation of rate coefficients.
In a recent paper, our group analyzed the influence of multiple conformers and paths in the evaluation of rate constants and relative abundances of products formed in the thermal decomposition of 1-propanol radicals using different methodologies including tsscds [66]. Specifically, the most relevant pathways reported in the literature [115,116,117,118,119,120,121] are obtained with tsscds, except for the barrierless dissociation leading to propene + OH, since the present version of tsscds cannot handle this type of reactions. Of significance, an important number of reactant and TS conformers, not described in the previous studies, are obtained with tsscds.
A conformational reaction channel (CRC) was defined in our study [66] as the group of all the paths that connect the conformers of a given reactant with the corresponding TS conformers. The influence of these conformers on the rate constants and branchings ratios was investigated in detail [66]. To study such influence, the output of tsscds (families of CRCs) was fed into a computer program to treat torsional anharmonicity named Q2DTOR (also developed in our group) [122]. The results obtained with tsscds and Q2DTOR were finally employed to calculate variational transition state theory (VTST) [123,124,125] rate constants for all the CRCs. The multipath (MP) approach within VTST was employed [125,126,127,128,129], where the rate constant of a given CRC is calculated using contributions from all the conformers and paths. For comparison purposes the simplest one-well (1W) approach is also considered; in the 1W method only the most stable conformers of reactant and TS are considered. As seen in Figure 7, the product abundances obtained in the temperature range 1000–2000 K are greatly influenced by the selected approach (MP vs 1W), particularly for the major products: ethene + CH2OH and formaldehyde + ethyl radical [66]. Our results show the importance of using automated codes for discovering reaction mechanisms and sampling potential energy surfaces.
Very recently, Fenard et al. developed a detailed kinetic model of the low-temperature oxidation of tetrahydrofuran (THF) based on theoretically-calculated rate constants [67]. The reaction pathways involved in these processes were probed with our automated software tsscds [67] using CBS-QB3 as the choice for the high-level of electronic structure. The rate constants were determined using TST with a tunneling correction using an Eckart potential.
The predictions from the model developed by Fenard et al. are overall in good agreement with the different experimental measurements. Namely, it reproduces ignition delay times obtained in a rapid-compression machine and in a shock tube, as well as numerous product mole fractions measured in a jet-stirred reactor.

3.4. Organometallic Catalysis

Computational studies of organometallic catalysis are becoming increasingly more important because they can help elucidate reaction mechanisms, characterize catalytic intermediates, supplement experimental studies, and also because of their predictive power [124,130,131,132,133].
However, the traditional workflow of most computational studies consists of using chemical intuition in the design of reaction routes and construction of guess TS structures. In recent years the appearance of powerful automated computational methods to study homogenous catalysis [27,43,134,135,136] very much eased the tedious work of manual searches.
To exemplify the use of tsscds in organometallic catalysis, the cobalt-catalyzed hydroformylation of ethylene was chosen [43]. Very briefly, the first step in our computational study was to generate all combinations of the catalyst Co(CO)3 with any of the starting materials (CO, H2 and ethylene), which in this case amounts to eight. Each of these combinations has fewer atoms than the overall system and they were named sub-systems in our original paper [43]. Standard tsscds is then run in each sub-system to build the reaction networks. Finally, the full reaction network is obtained after merging the individual results for each sub-system.
Figure 8 shows the tsscds-calculated free energy profile for the formation of propanal (C3H6O), which is the predominant channel; the level of theory employed was B3LYP/6-31G(d,p). As pointed out in the original paper, this is not the best electronic structure method for this system and it was only selected for comparison purposes. Additionally, we simulated the reactivity in the gas phase because, for this system, solvent effects are unimportant [43,133].
The mechanism shown in Figure 8 was obtained in an automated manner, and agrees with the one predict by Heck and Breslow in the 1960s [137] and with more recent mechanistic studies [133]. This is a very interesting result as we needed to make no assumptions in our automated calculations. Additionally, our method predicts that hydrogenation of ethylene is a side reaction that can be predominant under low CO partial pressures.
With the full reaction network constructed, the kinetics simulation module of tsscds can provide a rate law for the hydroformylation reaction when a range of different initial conditions for each species is employed. The kinetics calculations consist of transition state theory calculations [59,60,61,62] for the thermal rate constants at 423 K, and subsequent Monte Carlo simulations using different initial conditions of the reactants. Table 2 shows the orders of the catalyst and starting materials for the hydroformylation reaction obtained experimentally [138], with tsscds [43], using a kinetic model based on highly-accurate electronic structure calculations by Harvey and co-workers [133], and obtained from another automated method by Habershon [27].
As seen in Table 2, tsscds agrees rather well with experiment and with the results obtained by Harvey and co-workers [133]. Moreover, tsscds agrees much better with experiment than the other automated method does [27] (last column of Table 2), despite the fact that both employ the same alkene, initial conditions for the kinetics, and level of theory for the electronic structure calculations.

4. Improvements

In this section we describe some improvements we plan to implement in the near future. They include: the use of Spectral Graph Theory, implementation of knowledge-based methods, implementation of rare event acceleration MD simulations, interface with other electronic structure codes, reparametrization of semiempirical methods, and the study of condensed phase reactions.

4.1. Use of Spectral Graph Theory to Minimize the Number of Hessian Calculations

In standard tsscds, every single structure obtained after the BBFS analysis is subjected to TS optimization [45]. As seen in Figure 9a, for a trajectory i , BBFS selects m i TS candidates, which results in M = i = 1 n m i optimizations, where n is the total number of trajectories. On the one hand, these M optimizations are the most CPU-time consuming step of the procedure as they involve Hessian calculations, while the integration of the trajectories only requires gradients. On the other hand, a number of those optimizations are repeated. This is so because trajectories visit more often those areas of the configurational space around the kinetically most relevant TSs, leading to multiple optimizations of those structures.
The workflow of the enhanced procedure is shown in Figure 9b. Briefly, instead of carrying out the optimizations for every single structure selected by the BBFS algorithm (as in the original implementation), the new procedure will run the MD simulations and store at once the M structures for the analysis of all trajectory data. This analysis will consist of a pre-screening, a Spectral Graph Theory (SGT) step, and the final optimization step.
Upon completion of the MD simulations, a pre-screening of the M structures will be performed based on the eigenvalues of the Laplacian matrix [44]. As pointed out above, the lowest eigenvalues of this matrix indicate the degree of fragmentation of the molecular system. We aim here to discard highly fragmented structures, i.e., TSs connecting van der Waals complexes, usually of negligible relevance in a kinetics study. In the SGT step the remaining points will be partitioned into N groups according to the eigenvalues of a TS adjacency matrix, calculated as the average of the reactant and product adjacency matrices. Finally, we will select the closest point (geometry) to the centroid of each cluster for optimization. With this new scheme the gain in efficiency can easily be quantified as the reduction in the number of optimizations from M to N .

4.2. Implementation of Knowledge-Based Mechanism Generators

A number of reaction discovery methods are based on the so-called chemical heuristics [23,48,49,50]. In these methods, molecules are typically represented as graphs, in pretty much the same way as in tsscds. Then, by applying transformations, based on encoded rules or principles inspired by organic chemistry, to the reactant molecule graph, reactions, products and intermediates can readily be obtained. Compared to MD-based methods, heuristic-based methods are less CPU-time demanding.
Our idea will be to combine a heuristic-based bias in the MD simulations alongside with our BBFS algorithm to obtain TSs. In particular, having defined a set of encoded rules based on chemical knowledge, every single MD simulation will suffer a different bias, aimed to trigger a particular reaction mechanism. In this way, the problem of multiple optimizations of a given TS mentioned above would be minimized, if not completely avoided. The bias (analytical) potentials will be added on top of the semiempirical potential to steer the dynamics towards a particular intermediate or product.

4.3. Implementation of Rare-Event Acceleration MD Methods

One of the shortcomings of tsscds is the fact that chemical reactions are triggered by using very high energies in the MD simulations. While this approach was successfully employed to tackle different problems, it is biased towards the entropically favored reaction pathways. To alleviate this drawback of the method we propose to replace the current MD strategy by the rare-event acceleration method named Boxed Molecular Dynamics (BXD) [139]. BXD has its roots in work done by one of us and D. Shalashilin more than a decade ago [140]. It introduces several reflective barriers in the phase space of a MD trajectory along a particular collective variable. Those boundaries are employed to push the dynamics along the collective variable into regions of phase space which would be rarely sampled in an unbiased trajectory. However, the use of BXD constrains in configuration space suffers from the same “entropic” bias mentioned above.
A generalization of BXD has been very recently put forward by Glowacki and co-workers [141]. They show that the BXD bias can also be introduced along the potential energy (E) of the system, which is referred to as BXDE. By scanning through potential energy “boxes”, the energetic “windows” at which different chemical reaction channels switch on or off can be identified. The software design of tsscds is highly modular, which means that interfacing it with BXDE only requires little effort, like the need of compatible input/output geometry formats in both codes and the use of extra keywords in tsscds.

4.4. Interface with Other Electronic Structure Codes

At present tsscds has been only interfaced with the MOPAC2016 [75] and the G09 [76] electronic structure packages. The MD simulation employs gradients calculated at the semiempirical level of theory, and the optimization step is carried out at both the semiempirical level with MOPAC2016 and using higher levels (ab initio/DFT) with G09. Although we plan to reparametrize a semiempirical Hamiltonian for use in organometallic catalysis (see below), we do not want to be limited to this low-level electronic structure calculations. Therefore, we will use the ASE package [142] to interface tsscds with other electronic structure codes like NWCHEM [143] or ORCA [144].

4.5. Reparametrization of Semiempirical Methods

The application of the tsscds method relies on the use of semiempirical Hamiltonians for exploring potential energy surfaces. For this reason, it is important that the semiempirical method provides a reasonably accurate representation of the system under investigation. Although significant improvements in these methods have been made over the last years [145], there are still known limitations, which claim for further developments and more accurate parametrizations. Two important limitations concern the non-covalent interactions for large systems and ligand dissociation energies for transition metal complexes. In both cases, the performance of the semiempirical methods is, in general, quite poor. Our goal is therefore to improve the description of both non-covalent interactions and transition metal complexes in PM7.
Regarding non-covalent interactions, we aim to develop an analytical correction for PM7. To this end, we will consider a set of small molecules, which are representative of the most important functional groups. All pairs of molecules will be considered to calculate interaction energies at three levels of theory: coupled-cluster (CC), DFT and PM7. For every pair, various orientations will be considered, each one emphasizing a different two-body interaction.
Then, sums of two-body Buckingham potentials (supplemented with damping functions for the dispersion) will be fit to the CC, DFT and PM7 interaction energies using our genetic algorithm program GAFit [146]. Finally, the resulting potentials V fit , CC , V fit , DFT and V fit , PM 7 will be employed to build corrections V X c o r r to the PM7 interaction energies:
V X c o r r = V fit , X V fit , PM 7
where X is either CC or DFT. Whereas the V DFT c o r r correction term will be employed to validate this methodology as explained below, the highly-accurate V CC c o r r correction will be used once the validation succeeds.
The correction will be added to the PM7 energy V PM 7 so that the PM7 Hamiltonian corrected for non-covalent ( n c ) interactions would read:
V PM 7 , X n c = V PM 7 + V X c o r r
The strategy of using small representative molecules and sums of two-body functions was successfully employed in the development of intermolecular potentials for interactions of protonated peptides and silyl ions with perfluoroalkane self-assembled monolayers [147,148]. Nevertheless, this strategy will be validated for the new functional groups by running DFT calculations for large systems. This will allow us to compare the DFT-calculated energies with those obtained with V PM 7 , DFT n c .
The semiempirical methods, and particularly PM6 and PM7, do not perform well for transition-metal complexes [149]. Our strategy here will be to reoptimize the PM7 Hamiltonian as in previous studies of our group (e.g., see ref. [68]). We will select popular transition metals and ligand molecules used in organometallic catalysis, and will carry out high-level ab initio calculations for our own benchmark database. To gain flexibility in the parametrizations, we will consider the possibility of defining “atom types” for the ligand atoms, depending on the functional groups, in much the same way as that done for the parametrization of the hpCADD NDDO Hamiltonian [150].

4.6. Study of Condensed Phase Reactions

Our method is not limited to gas phase reactions. Although currently it only handles reactions in the gas phase, its modular design allows for a smooth adaptation of tsscds to deal with condensed phase reactions. For instance, to study solvent effects, the easiest way would be to use an implicit model, which in practice would only entail adding the appropriate keywords to the templates employed for the different electronic structure programs.
On the contrary, if one wants to use explicit solvent molecules, the MD module must be changed or substituted. At present, the MD module is a modified version of DRC routine in MOPAC2016, which includes different strategies for enhanced sampling, as detailed in the tutorial of tsscds [47]. To include solvent molecules in the MD simulations, one possibility would be to use CHARMM [151] or to adapt DRC. Finally, if the interest is a gas surface reactions, VENUS [152] would be the choice to run the MD simulations because the authors have vast experience using this program.

Author Contributions

Writing-Review & Editing, S.A.V., X.L.O. and E.M.-N.

Funding

This research was funded by “Consellería de Cultura, Educación e Ordenación Universitaria, Xunta de Galicia”, grant number ED431C 2017/17)”, and by “Ministerio de Economía y Competitividad of Spain”, grant number CTQ2014-58617-R.

Acknowledgments

The authors thank “Centro de Supercomputación de Galicia (CESGA)” for the use of their computational facilities.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

References

  1. Schlegel, H.B. Geometry optimization. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 790–809. [Google Scholar] [CrossRef]
  2. Davis, H.L.; Wales, D.J.; Berry, R.S. Exploring potential energy surfaces with transition state calculations. J. Chem. Phys. 1990, 92, 4308–4319. [Google Scholar] [CrossRef]
  3. Sun, J.Q.; Ruedenberg, K. Gradient Extremals and Steepest Descent Lines on Potential Energy Surfaces. J. Chem. Phys. 1993, 98, 9707–9714. [Google Scholar] [CrossRef]
  4. Tsai, C.J.; Jordan, K.D. Use of an eigenmode method to locate the stationary points on the potential energy surfaces of selected argon and water clusters. J. Phys. Chem. 1993, 97, 11227–11237. [Google Scholar] [CrossRef]
  5. Abashkin, Y.; Russo, N. Transition state structures and reaction profiles from constrained optimization procedure. Implementation in the framework of density functional theory. J. Chem. Phys. 1994, 100, 4477–4483. [Google Scholar] [CrossRef]
  6. Bondensgard, K.; Jensen, F. Gradient Extremal Bifurcation and Turning Points: An Application to the H2CO Potential Energy Surface. J. Chem. Phys. 1996, 104, 8025–8031. [Google Scholar] [CrossRef]
  7. Doye, J.P.K.; Wales, D.J. Surveying a potential energy surface by eigenvector-following. Z. Phys. D 1997, 40, 194–197. [Google Scholar] [CrossRef]
  8. Quapp, W.; Hirsch, M.; Imig, O.; Heidrich, D. Searching for Saddle Points of Potential Energy Surfaces by Following a Reduced Gradient. J. Comput. Chem. 1998, 19, 1087–1100. [Google Scholar] [CrossRef]
  9. Černohorský, M.; Kettou, S.; Koča, J. VADER:  New Software for Exploring Interconversions on Potential Energy Surfaces. J. Chem. Inf. Comput. Sci. 1999, 39, 705–712. [Google Scholar] [CrossRef]
  10. Westerberg, K.M.; Floudas, C.A. Locating all transition states and studying the reaction pathways of potential energy surfaces. J. Chem. Phys. 1999, 110, 9259–9295. [Google Scholar] [CrossRef]
  11. Wales, D.J.; Doye, J.P.; Miller, M.A.; Mortenson, P.N.; Walsh, T.R. Energy Landscapes: From Clusters to Biomolecules. Adv. Chem. Phys. 2000, 115, 1–111. [Google Scholar]
  12. Irikura, K.K.; Johnson, R.D. Predicting unexpected chemical reactions by isopotential searching. J. Phys. Chem. A 2000, 104, 2191–2194. [Google Scholar] [CrossRef]
  13. Müller, E.M.; Meijere, A.D.; Grubmüller, H. Predicting unimolecular chemical reactions: Chemical flooding. J. Chem. Phys. 2002, 116, 897–905. [Google Scholar] [CrossRef]
  14. Dallos, M.; Lischka, H.; Ventura Do Monte, E.; Hirsch, M.; Quapp, W. Determination of Energy Minima and Saddle Points Using Multireference Configuration Interaction Methods in Combination with Reduced Gradient Following: The S0 surface of H2CO and the T1 and T2 surfaces of acetylene. J. Comput. Chem. 2002, 23, 576–583. [Google Scholar] [CrossRef] [PubMed]
  15. Baker, J.; Wolinski, K. Isomerization of stilbene using enforced geometry optimization. J. Comput. Chem. 2011, 32, 43–53. [Google Scholar] [CrossRef] [PubMed]
  16. Zimmerman, P.M. Automated discovery of chemically reasonable elementary reaction steps. J. Comput. Chem. 2013, 34, 1385–1392. [Google Scholar] [CrossRef] [PubMed]
  17. Zimmerman, P.M. Growing string method with interpolation and optimization in internal coordinates: Method and examples. J. Chem. Phys. 2013, 138, 184102. [Google Scholar] [CrossRef] [PubMed]
  18. Zimmerman, P. Reliable Transition State Searches Integrated with the Growing String Method. J. Chem. Theory Comput. 2013, 9, 3043–3050. [Google Scholar] [CrossRef] [PubMed]
  19. Zimmerman, P.M. Single-ended transition state finding with the growing string method. J. Comput. Chem. 2015, 36, 601–611. [Google Scholar] [CrossRef] [PubMed]
  20. Zimmerman, P.M. Navigating molecular space for reaction mechanisms: An efficient, automated procedure. Mol. Simul. 2015, 41, 43–54. [Google Scholar] [CrossRef]
  21. Jafari, M.; Zimmerman, P.M. Reliable and efficient reaction path and transition state finding for surface reactions with the growing string method. J. Comput. Chem. 2017, 38, 645–658. [Google Scholar] [CrossRef] [PubMed]
  22. Dewyer, A.L.; Zimmerman, P.M. Finding reaction mechanisms, intuitive or otherwise. Org. Biomol. Chem. 2017, 15, 501–504. [Google Scholar] [CrossRef] [PubMed]
  23. Rappoport, D.; Galvin, C.J.; Zubarev, D.Y.; Aspuru-Guzik, A. Complex Chemical Reaction Networks from Heuristics-Aided Quantum Chemistry. J. Chem. Theory Comput. 2014, 10, 897–907. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Schaefer, B.; Mohr, S.; Amsler, M.; Goedecker, S. Minima hopping guided path search: An efficient method for finding complex chemical reaction pathways. J. Chem. Phys. 2014, 140, 214102. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Wales, D.J. Perspective: Insight into reaction coordinates and dynamics from the potential energy landscape. J. Chem. Phys. 2015, 142, 130901. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  26. Habershon, S. Sampling reactive pathways with random walks in chemical space: Applications to molecular dissociation and catalysis. J. Chem. Phys. 2015, 143, 094106. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Habershon, S. Automated Prediction of Catalytic Mechanism and Rate Law Using Graph-Based Reaction Path Sampling. J. Chem. Theory Comput. 2016, 12, 1786–1798. [Google Scholar] [CrossRef] [PubMed]
  28. Zhang, X.-J.; Liu, Z.-P. Reaction sampling and reactivity prediction using the stochastic surface walking method. Phys. Chem. Chem. Phys. 2015, 17, 2757–2769. [Google Scholar] [CrossRef] [PubMed]
  29. Wang, L.-P.; McGibbon, R.T.; Pande, V.S.; Martinez, T.J. Automated Discovery and Refinement of Reactive Molecular Dynamics Pathways. J. Chem. Theory Comput. 2016, 12, 638–649. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  30. Wang, L.-P.; Titov, A.; McGibbon, R.; Liu, F.; Pande, V.S.; Martínez, T.J. Discovering chemistry with an ab initio nanoreactor. Nat. Chem. 2014, 6, 1044–1048. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  31. Yang, M.; Zou, J.; Wang, G.; Li, S. Automatic Reaction Pathway Search via Combined Molecular Dynamics and Coordinate Driving Method. J. Phys. Chem. A 2017, 121, 1351–1361. [Google Scholar] [CrossRef] [PubMed]
  32. Jacobson, L.D.; Bochevarov, A.D.; Watson, M.A.; Hughes, T.F.; Rinaldo, D.; Ehrlich, S.; Steinbrecher, T.B.; Vaitheeswaran, S.; Philipp, D.M.; Halls, M.D.; et al. Automated Transition State Search and Its Application to Diverse Types of Organic Reactions. J. Chem. Theory Comput. 2017, 13, 5780–5797. [Google Scholar] [CrossRef] [PubMed]
  33. Ohno, K.; Maeda, S. A Scaled Hypersphere Search Method for the Topography of Reaction Pathways on the Potential Energy Surface. Chem. Phys. Lett. 2004, 384, 277–282. [Google Scholar] [CrossRef]
  34. Maeda, S.; Ohno, K. Global Mapping of Equilibrium and Transition Structures on Potential Energy Surfaces by the Scaled Hypersphere Search Method: Applications to ab Initio Surfaces of Formaldehyde and Propyne Molecules. J. Phys. Chem. A 2005, 109, 5742–5753. [Google Scholar] [CrossRef] [PubMed]
  35. Ohno, K.; Maeda, S. Global Reaction Route Mapping on Potential Energy Surfaces of Formaldehyde, Formic Acid, and Their Metail-Substituted Analogues. J. Phys. Chem. A 2006, 110, 8933–8941. [Google Scholar] [CrossRef] [PubMed]
  36. Ohno, K.; Maeda, S. Automated Exploration of Reaction Channels. Phys. Scr. 2008, 78, 058122. [Google Scholar] [CrossRef]
  37. Maeda, S.; Morokuma, K. Communications: A systematic method for locating transition structures of A + B → X type reactions. J. Chem. Phys. 2010, 132, 241102. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Maeda, S.; Morokuma, K. Finding Reaction Pathways of Type A + B → X: Toward Systematic Prediction of Reaction Mechanisms. J. Chem. Theory Comput. 2011, 7, 2335–2345. [Google Scholar] [CrossRef] [PubMed]
  39. Maeda, S.; Ohno, K.; Morokuma, K. Systematic exploration of the mechanism of chemical reactions: The global reaction route mapping (GRRM) strategy using the ADDF and AFIR methods. Phys. Chem. Chem. Phys. 2013, 15, 3683–3701. [Google Scholar] [CrossRef] [PubMed]
  40. Maeda, S.; Taketsugu, T.; Morokuma, K. Exploring transition state structures for intramolecular pathways by the artificial force induced reaction method. J. Comput. Chem. 2014, 35, 166–173. [Google Scholar] [CrossRef] [PubMed]
  41. Maeda, S.; Harabuchi, Y.; Takagi, M.; Taketsugu, T.; Morokuma, K. Artificial Force Induced Reaction (AFIR) Method for Exploring Quantum Chemical Potential Energy Surfaces. Chem. Rec. 2016, 16, 2232–2248. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Maeda, S.; Harabuchi, Y.; Takagi, M.; Saita, K.; Suzuki, K.; Ichino, T.; Sumiya, Y.; Sugiyama, K.; Ono, Y. Implementation and performance of the artificial force induced reaction method in the GRRM17 program. J. Comput. Chem. 2017, 39, 233–250. [Google Scholar] [CrossRef] [PubMed]
  43. Varela, J.A.; Vazquez, S.A.; Martinez-Nunez, E. An automated method to find reaction mechanisms and solve the kinetics in organometallic catalysis. Chem. Sci. 2017, 8, 3843–3851. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Martínez-Núñez, E. An automated transition state search using classical trajectories initialized at multiple minima. Phys. Chem. Chem. Phys. 2015, 17, 14912–14921. [Google Scholar] [CrossRef] [PubMed]
  45. Martínez-Núñez, E. An automated method to find transition states using chemical dynamics simulations. J. Comput. Chem. 2015, 36, 222–234. [Google Scholar] [CrossRef] [PubMed]
  46. Rodríguez, A.; Rodríguez-Fernández, R.; Vázquez, S.A.; Barnes, G.L.; Stewart, J.J.P.; Martínez-Núñez, E. tsscds2018: A code for automated discovery of chemical reaction mechanisms and solving the kinetics. J. Comput. Chem. 2018, 39, 1922–1930. [Google Scholar]
  47. Transition State Search Using Chemical Dynamics Simulations. Available online: http://forge.cesga.es/wiki/g/tsscds/HomePage (accessed on 25 October 2018).
  48. Broadbelt, L.J.; Stark, S.M.; Klein, M.T. Computer Generated Pyrolysis Modeling: On-the-Fly Generation of Species, Reactions, and Rates. Ind. Eng. Chem. Res. 1994, 33, 790–799. [Google Scholar] [CrossRef]
  49. Matheu, D.M.; Dean, A.M.; Grenda, J.M.; Green, W.H. Mechanism Generation with Integrated Pressure Dependence:  A New Model for Methane Pyrolysis. J. Phys. Chem. A 2003, 107, 8552–8565. [Google Scholar] [CrossRef]
  50. Gao, C.W.; Allen, J.W.; Green, W.H.; West, R.H. Reaction Mechanism Generator: Automatic construction of chemical kinetic mechanisms. Comput. Phys. Commun. 2016, 203, 212–225. [Google Scholar] [CrossRef] [Green Version]
  51. Bhoorasingh, P.L.; West, R.H. Transition state geometry prediction using molecular group contributions. Phys. Chem. Chem. Phys. 2015, 17, 32173–32182. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Bhoorasingh, P.L.; Slakman, B.L.; Seyedzadeh Khanshan, F.; Cain, J.Y.; West, R.H. Automated Transition State Theory Calculations for High-Throughput Kinetics. J. Phys. Chem. A 2017, 121, 6896–6904. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  53. Suleimanov, Y.V.; Green, W.H. Automated Discovery of Elementary Chemical Reaction Steps Using Freezing String and Berny Optimization Methods. J. Chem. Theory Comput. 2015, 11, 4248–4259. [Google Scholar] [CrossRef] [PubMed]
  54. Bergeler, M.; Simm, G.N.; Proppe, J.; Reiher, M. Heuristics-Guided Exploration of Reaction Mechanisms. J. Chem. Theory Comput. 2015, 11, 5712–5722. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  55. Proppe, J.; Husch, T.; Simm, G.N.; Reiher, M. Uncertainty quantification for quantum chemical models of complex reaction networks. Faraday Discuss. 2016, 195, 497–520. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  56. Simm, G.N.; Reiher, M. Context-Driven Exploration of Complex Chemical Reaction Networks. J. Chem. Theor. Comput. 2017, 13, 6108–6119. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  57. Simm, G.N.; Reiher, M. Error-Controlled Exploration of Chemical Reaction Networks with Gaussian Processes. J. Chem. Theor. Comput. 2018, 14, 5238–5248. [Google Scholar] [CrossRef] [PubMed]
  58. Dewyer, A.L.; Argüelles, A.J.; Zimmerman, P.M. Methods for exploring reaction space in molecular systems. WIREs Comput. Mol. Sci. 2018, 8, e1354. [Google Scholar] [CrossRef]
  59. Eyring, H. The Activated Complex in Chemical Reactions. J. Chem. Phys. 1935, 3, 107–115. [Google Scholar] [CrossRef]
  60. Wigner, E. The transition state method. Trans. Faraday Soc. 1938, 34, 29–41. [Google Scholar] [CrossRef]
  61. Keck, J.C. Variational Theory of Reaction Rates. Adv. Chem. Phys. 1967, 13, 85–121. [Google Scholar]
  62. Pechukas, P. Dynamics of Molecular Collisions; Plenum: New York, NY, USA, 1976. [Google Scholar]
  63. Baker, J. An algorithm for the location of transition states. J. Comput. Chem. 1986, 7, 385–395. [Google Scholar] [CrossRef]
  64. Fukui, K. The Path of Chemical Reactions-The IRC Approach. Acc. Chem. Res. 1981, 14, 363–368. [Google Scholar] [CrossRef]
  65. Gillespie, D.T. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 1976, 22, 403–434. [Google Scholar] [CrossRef]
  66. Ferro-Costas, D.; Martínez-Núñez, E.; Rodríguez-Otero, J.; Cabaleiro-Lago, E.; Estévez, C.M.; Fernández, B.; Fernández-Ramos, A.; Vázquez, S.A. Influence of Multiple Conformations and Paths on Rate Constants and Product Branching Ratios. Thermal Decomposition of 1-Propanol Radicals. J. Phys. Chem. A 2018, 122, 4790–4800. [Google Scholar] [CrossRef] [PubMed]
  67. Fenard, Y.; Gil, A.; Vanhove, G.; Carstensen, H.-H.; Van Geem, K.M.; Westmoreland, P.R.; Herbinet, O.; Battin-Leclerc, F. A model of tetrahydrofuran low-temperature oxidation based on theoretically calculated rate constants. Combust. Flame 2018, 191, 252–269. [Google Scholar] [CrossRef] [Green Version]
  68. Wilhelm, M.J.; Martínez-Núñez, E.; González-Vázquez, J.; Vázquez, S.A.; Smith, J.M.; Dai, H.-L. Is Photolytic Production a Viable Source of HCN and HNC in Astrophysical Environments? A Laboratory-based Feasibility Study of Methyl Cyanoformate. Astrophys. J. 2017, 849, 15. [Google Scholar] [CrossRef]
  69. Perez-Soto, R.; Vazquez, S.A.; Martinez-Nunez, E. Photodissociation of acryloyl chloride at 193 nm: Interpretation of the product energy distributions, and new elimination pathways. Phys. Chem. Chem. Phys. 2016, 18, 5019–5026. [Google Scholar] [CrossRef] [PubMed]
  70. Vazquez, S.A.; Martinez-Nunez, E. HCN elimination from vinyl cyanide: Product energy partitioning, the role of hydrogen-deuterium exchange reactions and a new pathway. Phys. Chem. Chem. Phys. 2015, 17, 6948–6955. [Google Scholar] [CrossRef] [PubMed]
  71. Rossich Molina, E.; Salpin, J.-Y.; Spezia, R.; Martinez-Nunez, E. On the gas phase fragmentation of protonated uracil: A statistical perspective. Phys. Chem. Chem. Phys. 2016, 18, 14980–14990. [Google Scholar] [CrossRef] [PubMed]
  72. Simm, G.N.; Vaucher, A.C.; Reiher, M. Exploration of Reaction Pathways and Chemical Transformation Networks. J. Phys. Chem. A 2018. [Google Scholar] [CrossRef] [PubMed]
  73. Hase, W.L.; Buckowski, D.G. Monte carlo sampling of a microcanonical ensemble of classical harmonic oscillators. Chem. Phys. Lett. 1980, 74, 284–287. [Google Scholar] [CrossRef]
  74. Bougueroua, S.; Spezia, R.; Pezzotti, S.; Vial, S.; Quessette, F.; Barth, D.; Gaigeot, M.-P. Graph theory for automatic structural recognition in molecular dynamics simulations. J. Chem. Phys. 2018, 149, 184102. [Google Scholar] [CrossRef] [PubMed]
  75. Stewart, J.J.P. MOPAC2016, Stewart Computational Chemistry. Available online: http://openmopac.net (accessed on 20 October 2018).
  76. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G.A.; et al. Gaussian 09 Revision A.02; Gaussian Inc.: Wallingford, CT, USA, 2009. [Google Scholar]
  77. Pietrucci, F.; Andreoni, W. Graph Theory Meets Ab Initio Molecule Dynamics: Atomic Structures and Transformations at the Nanoscale. Phys. Rev. Lett. 2011, 107, 085504. [Google Scholar] [CrossRef] [PubMed]
  78. Smith, G.; Gilbert, R.G. Theory of Unimolecular and Recombination Reactions; Blackwell Scientific Publications: Oxford, UK, 1990. [Google Scholar]
  79. Tarrazo-Antelo, T.; Martinez-Nunez, E.; Vazquez, S.A. Ab initio and RRKM study of the elimination of HF and HCl from chlorofluoroethylene. Chem. Phys. Lett. 2007, 435, 176–181. [Google Scholar] [CrossRef]
  80. Martínez-Núñez, E.; Vázquez, S. Rotational distributions of HBr in the photodissociation of vinyl bromide at 193 nm: An investigation by direct quasiclassical trajectory calculations. Chem. Phys. Lett. 2006, 425, 22–27. [Google Scholar] [CrossRef]
  81. Martínez-Núñez, E.; Vázquez, S. Quasiclassical trajectory calculations on the photodissociation of CF2CHCl at 193 nm: Product energy distributions for the HF and HCl eliminations. J. Chem. Phys. 2005, 122, 104316. [Google Scholar] [CrossRef] [PubMed]
  82. Martínez-Núñez, E.; Vázquez, S.A.; Aoiz, F.J.; Bañares, L.; Castillo, J.F. Further investigation of the HCl elimination in the photodissociation of vinyl chloride at 193 nm: A direct MP2/6-31G(d,p) trajectory study. Chem. Phys. Lett. 2004, 386, 225–232. [Google Scholar] [CrossRef]
  83. Martínez-Núñez, E.; Vázquez, S. Rovibrational distributions of HF in the photodissociation of vinyl fluoride at 193 nm: A direct MP2 quasiclassical trajectory study. J. Chem. Phys. 2004, 121, 5179–5182. [Google Scholar] [CrossRef] [PubMed]
  84. Martínez-Núñez, E.; Fernández-Ramos, A.; Vázquez, S.A.; JavierAoiz, F.; Bañares, L. A Direct Classical Trajectory Study of HCl Elimination from the 193 nm Photodissociation of Vinyl Chloride. J. Phys. Chem. A 2003, 107, 7611–7618. [Google Scholar] [CrossRef]
  85. Gonzalez-Vazquez, J.; Martinez-Nunez, E.; Fernandez-Ramos, A.; Vazquez, S.A. Dissociation of difluoroethylenes. II. Direct Classical Trajectory Study of the HF elimination from 1,2-difluoroethylene. J. Phys. Chem. A 2003, 107, 1398–1404. [Google Scholar] [CrossRef]
  86. Gonzalez-Vazquez, J.; Fernandez-Ramos, A.; Martinez-Nunez, E.; Vazquez, S.A. Dissociation of difluoroethylenes. I. Global potential energy surface, RRKM, and VTST calculations. J. Phys. Chem. A 2003, 107, 1389–1397. [Google Scholar] [CrossRef]
  87. Martínez-Núñez, E.; Estévez, C.M.; Flores, J.R.; Vázquez, S.A. Product energy distributions for the four-center HF elimination from 1,1-difluoroethylene. a direct dynamics study. Chem. Phys. Lett. 2001, 348, 81–88. [Google Scholar] [CrossRef]
  88. Martínez-Núñez, E.; Vázquez, S.A. Three-center vs. four-center HF elimination from vinyl fluoride: A direct dynamics study. Chem. Phys. Lett. 2000, 332, 583–590. [Google Scholar] [CrossRef]
  89. Homayoon, Z.; Vázquez, S.A.; Rodríguez-Fernández, R.; Martínez-Núñez, E. Ab initio and RRKM study of the HCN/HNC elimination channels from vinyl cyanide. J. Phys. Chem. A 2011, 115, 979–985. [Google Scholar] [CrossRef] [PubMed]
  90. Martinez-Nunez, E.; Vazquez, S.A.; Borges, I.; Rocha, A.B.; Estevez, C.M.; Castillo, J.F.; Aoiz, F.J. On the conformational memory in the photodissociation of formic acid. J. Phys. Chem. A 2005, 109, 2836–2839. [Google Scholar] [CrossRef] [PubMed]
  91. Martinez-Nunez, E.; Vazquez, S.; Granucci, G.; Persico, M.; Estevez, C.M. Photodissociation of formic acid: A trajectory surface hopping study. Chem. Phys. Lett. 2005, 412, 35–40. [Google Scholar] [CrossRef]
  92. Chang, C.M.; Huang, Y.H.; Liu, S.Y.; Lee, Y.P.; Pombar-Perez, M.; Martinez-Nunez, E.; Vazquez, S.A. Internal energy of HCl upon photolysis of 2-chloropropene at 193 nm investigated with time-resolved Fourier-transform spectroscopy and quasiclassical trajectories. J. Chem. Phys. 2008, 129, 224301. [Google Scholar] [CrossRef] [PubMed]
  93. Spezia, R.; Martínez-Nuñez, E.; Vazquez, S.; Hase, W.L. Theoretical and computational studies of non-equilibrium and non-statistical dynamics in the gas phase, in the condensed phase and at interfaces. Philos. Trans. R. Soc. A 2017, 375, 20170035. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  94. Tsutsumi, T.; Harabuchi, Y.; Ono, Y.; Maeda, S.; Taketsugu, T. Analyses of trajectory on-the-fly based on the global reaction route map. Phys. Chem. Chem. Phys. 2018, 20, 1364–1372. [Google Scholar] [CrossRef] [PubMed]
  95. Wilhelm, M.J.; Nikow, M.; Letendre, L.; Dai, H.-L. Photodissociation of vinyl cyanide at 193 nm: Nascent product distributions of the molecular elimination channels. J. Chem. Phys. 2009, 130, 044307. [Google Scholar] [CrossRef] [PubMed]
  96. Chin, C.-H.; Lee, S.-H. Theoretical study of isomerization and decomposition of propenal. J. Chem. Phys. 2011, 134, 044309. [Google Scholar] [CrossRef] [PubMed]
  97. Chaudhuri, C.; Lee, S.-H. A complete look at the multi-channel dissociation of propenal photoexcited at 193 nm: Branching ratios and distributions of kinetic energy. Phys. Chem. Chem. Phys. 2011, 13, 7312–7321. [Google Scholar] [CrossRef] [PubMed]
  98. Lee, P.-W.; Scrape, P.G.; Butler, L.J.; Lee, Y.-P. Two HCl-Elimination Channels and Two CO-Formation Channels Detected with Time-Resolved Infrared Emission upon Photolysis of Acryloyl Chloride [CH2CHC(O)Cl] at 193 nm. J. Phys. Chem. A 2015, 119, 7293–7304. [Google Scholar] [CrossRef] [PubMed]
  99. Bauer, C.A.; Grimme, S. How to Compute Electron Ionization Mass Spectra from First Principles. J. Phys. Chem. A 2016, 120, 3755–3766. [Google Scholar] [CrossRef] [PubMed]
  100. Macaluso, V.; Homayoon, Z.; Spezia, R.; Hase, W.L. Threshold for shattering fragmentation in collision-induced dissociation of the doubly protonated tripeptide TIK(H+)2. Phys. Chem. Chem. Phys. 2018, 20, 19744–19749. [Google Scholar] [CrossRef] [PubMed]
  101. Martin-Somer, A.; Martens, J.; Grzetic, J.; Hase, W.L.; Oomens, J.; Spezia, R. Unimolecular Fragmentation of Deprotonated Diproline [Pro2-H] Studied by Chemical Dynamics Simulations and IRMPD Spectroscopy. J. Phys. Chem. A 2018, 122, 2612–2625. [Google Scholar] [CrossRef] [PubMed]
  102. Homayoon, Z.; Macaluso, V.; Martin-Somer, A.; Muniz, M.C.N.B.; Borges, I.; Hase, W.L.; Spezia, R. Chemical dynamics simulations of CID of peptide ions: Comparisons between TIK(H+)2 and TLK(H+)2 fragmentation dynamics, and with thermal simulations. Phys. Chem. Chem. Phys. 2018, 20, 3614–3629. [Google Scholar] [CrossRef] [PubMed]
  103. Martin-Somer, A.; Spezia, R.; Yáñez, M. Gas-phase reactivity of [Ca(formamide)]2+ complex: An example of different dynamical behaviours. Philos. Trans. R. Soc. A 2017, 375, 20160196. [Google Scholar] [CrossRef] [PubMed]
  104. Molina, E.R.; Eizaguirre, A.; Haldys, V.; Urban, D.; Doisneau, G.; Bourdreux, Y.; Beau, J.-M.; Salpin, J.-Y.; Spezia, R. Characterization of Protonated Model Disaccharides from Tandem Mass Spectrometry and Chemical Dynamics Simulations. ChemPhysChem 2017, 18, 2812–2823. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  105. Lee, G.; Park, E.; Chung, H.; Jeanvoine, Y.; Song, K.; Spezia, R. Gas phase fragmentation mechanisms of protonated testosterone as revealed by chemical dynamics simulations. Int. J. Mass Spectrom. 2016, 407, 40–50. [Google Scholar] [CrossRef]
  106. Spezia, R.; Lee, S.B.; Cho, A.; Song, K. Collision-induced dissociation mechanisms of protonated penta- and octa-glycine as revealed by chemical dynamics simulations. Int. J. Mass Spectrom. 2015, 392, 125–138. [Google Scholar] [CrossRef]
  107. Spezia, R.; Martens, J.; Oomens, J.; Song, K. Collision-induced dissociation pathways of protonated Gly2NH2 and Gly3NH2 in the short time-scale limit by chemical dynamics and ion spectroscopy. Int. J. Mass Spectrom. 2015, 388, 40–52. [Google Scholar] [CrossRef]
  108. Song, K.; Spezia, R. Theoretical Mass Spectrometry, Tracing Ions with Classical Trajectories; De Gruyter: Berlin, Germany, 2018. [Google Scholar]
  109. Pratihar, S.; Barnes, G.L.; Laskin, J.; Hase, W.L. Dynamics of Protonated Peptide Ion Collisions with Organic Surfaces: Consonance of Simulation and Experiment. J. Phys. Chem. Lett. 2016, 7, 3142–3150. [Google Scholar] [CrossRef] [PubMed]
  110. Pratihar, S.; Barnes, G.L.; Hase, W.L. Chemical dynamics simulations of energy transfer, surface-induced dissociation, soft-landing, and reactive-landing in collisions of protonated peptide ions with organic surfaces. Chem. Soc. Rev. 2016, 45, 3595–3608. [Google Scholar] [CrossRef] [PubMed]
  111. Barnes, G.L.; Young, K.; Yang, L.; Hase, W.L. Fragmentation and reactivity in collisions of protonated diglycine with chemically modified perfluorinated alkylthiolate-self-assembled monolayer surfaces. J. Chem. Phys. 2011, 134, 094106. [Google Scholar] [CrossRef] [PubMed]
  112. Park, K.; Deb, B.; Song, K.; Hase, W.L. Importance of Shattering Fragmentation in the Surface-Induced Dissociation of Protonated Octaglycine. JASMS 2009, 20, 939–948. [Google Scholar] [CrossRef] [PubMed]
  113. Barnes, G.L.; Hase, W.L. Energy Transfer, Unfolding, and Fragmentation Dynamics in Collisions of N-Protonated Octaglycine with an H-SAM Surface. J. Am. Chem. Soc. 2009, 131, 17185–17193. [Google Scholar] [CrossRef] [PubMed]
  114. Martínez-Núñez, E.; Fernández-Ramos, A.; Vázquez, S.A.; Marques, J.M.C.; Xue, M.; Hase, W.L. Quasiclassical dynamics simulation of the collision-induced dissociation of Cr (CO)6 + with Xe. J. Chem. Phys. 2005, 123, 154311. [Google Scholar] [CrossRef] [PubMed]
  115. Zador, J.; Jasper, A.W.; Miller, J.A. The reaction between propene and hydroxyl. Phys. Chem. Chem. Phys. 2009, 11, 11040–11053. [Google Scholar] [CrossRef] [PubMed]
  116. Zhou, C.-W.; Li, Z.-R.; Li, X.-Y. Kinetics and Mechanism for Formation of Enols in Reaction of Hydroxide Radical with Propene. J. Phys. Chem. A 2009, 113, 2372–2382. [Google Scholar] [CrossRef] [PubMed]
  117. Huynh, L.K.; Zhang, H.R.; Zhang, S.; Eddings, E.; Sarofim, A.; Law, M.E.; Westmoreland, P.R.; Truong, T.N. Kinetics of Enol Formation from Reaction of OH with Propene. J. Phys. Chem. A 2009, 113, 3177–3185. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  118. El-Nahas, A.M.; Uchimaru, T.; Sugie, M.; Tokuhashi, K.; Sekiya, A. Relative reactivity and regioselectivity of halogen-substituted ethenes and propene toward addition of an OH radical or O (3P) atom: An ab initio study. THEOCHEM 2006, 770, 59–65. [Google Scholar] [CrossRef]
  119. Szori, M.; Fittschen, C.; Csizmadia, I.G.; Viskolcz, B. Allylic H-Abstraction Mechanism:  The Potential Energy Surface of the Reaction of Propene with OH Radical. J. Chem. Theor. Comput. 2006, 2, 1575–1586. [Google Scholar] [CrossRef] [PubMed]
  120. Díaz-Acosta, I.; Alvarez-Idaboy, J.R.; Vivier-Bunge, A. Mechanism of the OH-propene-O2 reaction: An ab initio study. Int. J. Chem. Kinet. 1999, 31, 29–36. [Google Scholar] [CrossRef]
  121. Alvarez-Idaboy, J.R.; Díaz-Acosta, I.; Vivier-Bunge, A. Energetics of mechanism of OH-propene reaction at low pressures in inert atmosphere. J. Comput. Chem. 1998, 19, 811–819. [Google Scholar] [CrossRef]
  122. Ferro-Costas, D.; Cordeiro, M.N.D.S.; Truhlar, D.G.; Fernández-Ramos, A. Q2DTor: A program to treat torsional anharmonicity through coupled pair torsions in flexible molecules. Comput. Phys. Commun. 2018, 232, 190–205. [Google Scholar] [CrossRef]
  123. Truhlar, D.G.; Isaacson, A.D.; Garret, G.C. Theory of Chemical Reaction Dynamics; Baer, M., Ed.; CRC: Boca Raton, FL, USA, 1985; Volume 4, p. 65. [Google Scholar]
  124. Schwarz, H. Chemistry with Methane: Concepts Rather than Recipes. Angew. Chem. Int. Ed. 2011, 50, 10096–10115. [Google Scholar] [CrossRef] [PubMed]
  125. Bao, J.L.; Truhlar, D.G. Variational transition state theory: Theoretical framework and recent developments. Chem. Soc. Rev. 2017, 46, 7548–7596. [Google Scholar] [CrossRef] [PubMed]
  126. Yu, T.; Zheng, J.; Truhlar, D.G. Multi-structural variational transition state theory. Kinetics of the 1,4-hydrogen shift isomerization of the pentyl radical with torsional anharmonicity. Chem. Sci. 2011, 2, 2199–2213. [Google Scholar] [CrossRef]
  127. Bao, J.L.; Meana-Pañeda, R.; Truhlar, D.G. Multi-path variational transition state theory for chiral molecules: The site-dependent kinetics for abstraction of hydrogen from 2-butanol by hydroperoxyl radical, analysis of hydrogen bonding in the transition state, and dramatic temperature dependence of the activation energy. Chem. Sci. 2015, 6, 5866–5881. [Google Scholar] [PubMed]
  128. Yu, T.; Zheng, J.; Truhlar, D.G. Multipath Variational Transition State Theory: Rate Constant of the 1,4-Hydrogen Shift Isomerization of the 2-Cyclohexylethyl Radical. J. Phys. Chem. A 2012, 116, 297–308. [Google Scholar] [CrossRef] [PubMed]
  129. Meana-Pañeda, R.; Fernández-Ramos, A. Accounting for conformational flexibility and torsional anharmonicity in the H + CH3CH2OH hydrogen abstraction reactions: A multi-path variational transition state theory study. J. Chem. Phys. 2014, 140, 174303. [Google Scholar] [CrossRef] [PubMed]
  130. Sperger, T.; Sanhueza, I.A.; Schoenebeck, F. Computation and Experiment: A Powerful Combination to Understand and Predict Reactivities. Acc. Chem. Res. 2016, 49, 1311–1319. [Google Scholar] [CrossRef] [PubMed]
  131. Peng, Q.; Paton, R.S. Catalytic Control in Cyclizations: From Computational Mechanistic Understanding to Selectivity Prediction. Acc. Chem. Res. 2016, 49, 1042–1051. [Google Scholar] [CrossRef] [PubMed]
  132. Sperger, T.; Sanhueza, I.A.; Kalvet, I.; Schoenebeck, F. Computational Studies of Synthetically Relevant Homogeneous Organometallic Catalysis Involving Ni, Pd, Ir, and Rh: An Overview of Commonly Employed DFT Methods and Mechanistic Insights. Chem. Rev. 2015, 115, 9532–9586. [Google Scholar] [CrossRef] [PubMed]
  133. Rush, L.E.; Pringle, P.G.; Harvey, J.N. Computational Kinetics of Cobalt-Catalyzed Alkene Hydroformylation. Angew. Chem. Int. Ed. 2014, 53, 8672–8676. [Google Scholar] [CrossRef] [PubMed]
  134. Maeda, S.; Morokuma, K. Toward Predicting Full Catalytic Cycle Using Automatic Reaction Path Search Method: A Case Study on HCo(CO)3-Catalyzed Hydroformylation. J. Chem. Theor. Comput. 2012, 8, 380–385. [Google Scholar] [CrossRef] [PubMed]
  135. Kim, Y.; Choi, S.; Kim, W.Y. Efficient Basin-Hopping Sampling of Reaction Intermediates through Molecular Fragmentation and Graph Theory. J. Chem. Theory Comput. 2014, 10, 2419–2426. [Google Scholar] [CrossRef] [PubMed]
  136. Kim, Y.; Kim, J.W.; Kim, Z.; Kim, W.Y. Efficient prediction of reaction paths through molecular graph and reaction network analysis. Chem. Sci. 2018, 9, 825–835. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  137. Heck, R.F.; Breslow, D.S. The Reaction of Cobalt Hydrotetracarbonyl with Olefins. J. Am. Chem. Soc. 1961, 83, 4023–4027. [Google Scholar] [CrossRef]
  138. Gholap, R.V.; Kut, O.M.; Bourne, J.R. Hydroformylation of propylene using an unmodified cobalt carbonyl catalyst: A kinetic study. Ind. Eng. Chem. Res. 1992, 31, 1597–1601. [Google Scholar] [CrossRef]
  139. Booth, J.; Vazquez, S.; Martínez-Núñez, E.; Marks, A.; Rodgers, J.; Glowacki, D.R.; Shalashilin, D.V. Recent Applications of Boxed Molecular Dynamics: A Simple Multiscale Technique for Atomistic Simulations. Philos. Trans. R. Soc. A 2014, 372, 20130384. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  140. Martínez-Núñez, E.; Shalashilin, D.V. Acceleration of classical mechanics by phase space constraints. J. Chem. Theor. Comput. 2006, 2, 912. [Google Scholar] [CrossRef] [PubMed]
  141. Shannon, R.J.; Amabilino, S.; O’Connor, M.; Shalishilin, D.V.; Glowacki, D.R. Adaptively Accelerating Reactive Molecular Dynamics Using Boxed Molecular Dynamics in Energy Space. J. Chem. Theor. Comput. 2018, 14, 4541–4552. [Google Scholar] [CrossRef] [PubMed]
  142. Larsen, A.H.; Mortensen, J.J.; Blomqvist, J.; Castelli, I.E.; Christensen, R.; Dułak, M.; Friis, J.; Groves, M.N.; Hammer, B.; Hargus, C.; et al. The atomic simulation environment—A Python library for working with atoms. J. Phys. Condens. Matter 2017, 29, 273002. [Google Scholar] [CrossRef] [PubMed]
  143. Valiev, M.; Bylaska, E.J.; Govind, N.; Kowalski, K.; Straatsma, T.P.; Van Dam, H.J.J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T.L.; et al. NWChem: A comprehensive and scalable open-source solution for large scale molecular simulations. Comput. Phys. Commun. 2010, 181, 1477–1489. [Google Scholar] [CrossRef] [Green Version]
  144. Neese, F. The ORCA program system. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 73–78. [Google Scholar] [CrossRef]
  145. Christensen, A.S.; Kubař, T.; Cui, Q.; Elstner, M. Semiempirical Quantum Mechanical Methods for Noncovalent Interactions for Chemical and Biochemical Applications. Chem. Rev. 2016, 116, 5301–5337. [Google Scholar] [CrossRef] [PubMed]
  146. Rodríguez-Fernández, R.; Pereira, F.B.; Marques, J.M.C.; Martínez-Núñez, E.; Vázquez, S.A. GAFit: A general-purpose, user-friendly program for fitting potential energy surfaces. Comput. Phys. Commun. 2017, 217, 89–98. [Google Scholar] [CrossRef]
  147. Nogueira, J.J.; Sánchez-Coronilla, A.; Marques, J.M.C.; Hase, W.L.; Martínez-Núñez, E.; Vázquez, S.A. Intermolecular potentials for simulations of collisions of SiNCS+ and (CH3)2SiNCS+ ions with fluorinated self-assembled monolayers. Chem. Phys. 2012, 399, 193–204. [Google Scholar] [CrossRef]
  148. Pratihar, S.; Kohale, S.C.; Vázquez, S.A.; Hase, W.L. Intermolecular Potential for Binding of Protonated Peptide Ions with Perfluorinated Hydrocarbon Surfaces. J. Phys. Chem. B 2014, 118, 5577–5588. [Google Scholar] [CrossRef] [PubMed]
  149. Semiempirical Molecular Orbital Models Based on the Neglect of Diatomic Differential Overlap Approximation. Available online: https://arxiv.org/abs/1806.06147 (accessed on 20 October 2018).
  150. Thomas, H.B.; Hennemann, M.; Kibies, P.; Hoffgaard, F.; Güssregen, S.; Hessler, G.; Kast, S.M.; Clark, T. The hpCADD NDDO Hamiltonian: Parametrization. J. Chem. Inf. Model. 2017, 57, 1907–1922. [Google Scholar] [CrossRef] [PubMed]
  151. Brooks, B.R.; Brooks, C.L.; Mackerell, A.D., Jr.; Nilsson, L.; Petrella, R.J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, S.; et al. CHARMM: The biomolecular simulation program. J. Comput. Chem. 2009, 30, 1545–1614. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  152. Hase, W.L.; Bolton, K.; Sainte Claire, P.D.; Duchovic, R.J.; Hu, X.; Komornicki, A.; Li, G.; Lim, K.F.; Lu, D.-H.; Peslherbe, G.H.; et al. Venus05: A General Chemical Dynamics Computer Program; Bloomington: Bloomington, IN, USA, 2004. [Google Scholar]
Figure 1. New HCN elimination mechanism from VCN obtained with tsscds. The numbers are relative energies (including the zero-point vibrational energy) with respect to VCN, calculated at the CCSD(T)/6-311++G(3df,3pd)//CCSD/6-311+G(2d,2p) level of theory with the vibrational frequencies obtained using CCSD/6-311+G(2d,2p) numerical Hessians.
Figure 1. New HCN elimination mechanism from VCN obtained with tsscds. The numbers are relative energies (including the zero-point vibrational energy) with respect to VCN, calculated at the CCSD(T)/6-311++G(3df,3pd)//CCSD/6-311+G(2d,2p) level of theory with the vibrational frequencies obtained using CCSD/6-311+G(2d,2p) numerical Hessians.
Molecules 23 03156 g001
Figure 2. Kinetic simulation results of the different HCN elimination channels from VCN.
Figure 2. Kinetic simulation results of the different HCN elimination channels from VCN.
Molecules 23 03156 g002
Figure 3. Minima obtained by tsscds for the C3H4O system. The structures are arranged in ascending order of their relative energies (shown at the bottom of each structure), which are obtained at the CCSD(T)/6-311+G(3df,2p)//B3LYP/6-311G(d,p) level of theory. Conformers are not included in the figure and only the lowest lying of each family is displayed.
Figure 3. Minima obtained by tsscds for the C3H4O system. The structures are arranged in ascending order of their relative energies (shown at the bottom of each structure), which are obtained at the CCSD(T)/6-311+G(3df,2p)//B3LYP/6-311G(d,p) level of theory. Conformers are not included in the figure and only the lowest lying of each family is displayed.
Molecules 23 03156 g003
Figure 4. Structure of AC minimum and the three new TSs found with tsscds for the HCl elimination from AC. Numbers are relative energies in kcal/mol (including the zero-point vibrational energy) with respect to AC, calculated at the CCSD(T)/6-311+G(3df,2p)//B3LYP/6-311+G(2d,2p) level of theory.
Figure 4. Structure of AC minimum and the three new TSs found with tsscds for the HCl elimination from AC. Numbers are relative energies in kcal/mol (including the zero-point vibrational energy) with respect to AC, calculated at the CCSD(T)/6-311+G(3df,2p)//B3LYP/6-311+G(2d,2p) level of theory.
Molecules 23 03156 g004
Figure 5. Relevant HCN and HNC pathways in the ground-state PES of methyl cyanoformate for an excitation energy of 148 kcal/mol. Relative energies (in kcal·mol−1) include ZPE contributions and were obtained by CCSD(T)/6-311++G(3df,3pd)//MP2/6-311+G(2d,2p) calculations.
Figure 5. Relevant HCN and HNC pathways in the ground-state PES of methyl cyanoformate for an excitation energy of 148 kcal/mol. Relative energies (in kcal·mol−1) include ZPE contributions and were obtained by CCSD(T)/6-311++G(3df,3pd)//MP2/6-311+G(2d,2p) calculations.
Molecules 23 03156 g005
Figure 6. Experimental (exp) and calculated (comp) intensities of precursor and fragment ions produced in the fragmentation of protonated uracil.
Figure 6. Experimental (exp) and calculated (comp) intensities of precursor and fragment ions produced in the fragmentation of protonated uracil.
Molecules 23 03156 g006
Figure 7. Branching ratios obtained in the kinetics simulations starting from one of the isomers of 1-propanol (only the two major mechanisms are shown). The solid and dashed lines correspond to the MP and 1W results, respectively.
Figure 7. Branching ratios obtained in the kinetics simulations starting from one of the isomers of 1-propanol (only the two major mechanisms are shown). The solid and dashed lines correspond to the MP and 1W results, respectively.
Molecules 23 03156 g007
Figure 8. Free energy profile for the Co-catalyzed hydroformylation of ethylene obtained in our tsscds study using DFT calculations [133].
Figure 8. Free energy profile for the Co-catalyzed hydroformylation of ethylene obtained in our tsscds study using DFT calculations [133].
Molecules 23 03156 g008
Figure 9. (a) Original tsscds showcasing an example with n different trajectories resulting in a total number of M = i = 1 n m i optimizations. (b) Modified tsscds showcasing the same example as in panel (a) with n different trajectories resulting in a total number of N optimizations.
Figure 9. (a) Original tsscds showcasing an example with n different trajectories resulting in a total number of M = i = 1 n m i optimizations. (b) Modified tsscds showcasing the same example as in panel (a) with n different trajectories resulting in a total number of N optimizations.
Molecules 23 03156 g009
Table 1. Relative product abundances obtained by different computational studies and experiment in the photodissociation of ACRL at 193 nm.
Table 1. Relative product abundances obtained by different computational studies and experiment in the photodissociation of ACRL at 193 nm.
ChannelChin et al. [96]TsscdsExp [97]
H2O0.010.030.07
CH2O0.650.200.07
H20.090.190.00
CO1.001.001.00
H2 + CO + HCCH6.821.491.10
Table 2. Orders of the hydroformylation reaction with respect to the catalyst and starting materials.
Table 2. Orders of the hydroformylation reaction with respect to the catalyst and starting materials.
SpeciesExp [138]tsscds [43]Harvey [133]Habershon [27]
H20.60.40.51
CO<0<0<0<0
catalyst0.80.50.51
alkene1110.55

Share and Cite

MDPI and ACS Style

Vázquez, S.A.; Otero, X.L.; Martinez-Nunez, E. A Trajectory-Based Method to Explore Reaction Mechanisms. Molecules 2018, 23, 3156. https://doi.org/10.3390/molecules23123156

AMA Style

Vázquez SA, Otero XL, Martinez-Nunez E. A Trajectory-Based Method to Explore Reaction Mechanisms. Molecules. 2018; 23(12):3156. https://doi.org/10.3390/molecules23123156

Chicago/Turabian Style

Vázquez, Saulo A., Xose L. Otero, and Emilio Martinez-Nunez. 2018. "A Trajectory-Based Method to Explore Reaction Mechanisms" Molecules 23, no. 12: 3156. https://doi.org/10.3390/molecules23123156

APA Style

Vázquez, S. A., Otero, X. L., & Martinez-Nunez, E. (2018). A Trajectory-Based Method to Explore Reaction Mechanisms. Molecules, 23(12), 3156. https://doi.org/10.3390/molecules23123156

Article Metrics

Back to TopTop