Advances in Automated Transition State Theory Calculations: Improvements on the AutoTST Framework

Kinetic modeling of combustion chemistry has made substantial progress in recent years with the development of increasingly detailed models. However, many of the chemical kinetic parameters utilized in detailed models are estimated, often inaccurately. To help replace rate estimates with more accurate calculations, we have developed AutoTST, an automated Transition State Theory rate calculator. This work describes improvements to AutoTST, including: a systematic conformer search to find an ensemble of low energy conformers, vibrational analysis to validate transition state geometries, more accurate symmetry number calculations, and a hindered rotor treatment when deriving kinetics. These improvements resulted in location of transition state geometry for 93% of cases and generation of kinetic parameters for 74% of cases. Newly calculated parameters agree well with benchmark calculations and perform well when used to replace estimated parameters in a detailed kinetic model of methanol combustion. Abstract Kinetic modeling of combustion chemistry has made substantial progress in recent years with the development of increasingly detailed models. However, many of the chemical kinetic parameters utilized in detailed models are estimated, often inaccurately. To help replace rate estimates with more accurate calculations, we have developed AutoTST, an automated Transition State Theory rate calculator. This work describes improvements to AutoTST, including: a systematic conformer search to ﬁnd an ensemble of low energy conformers, vibrational analysis to validate transition state geometries, more accurate symmetry number calculations, and a hindered rotor treatment when deriving kinetics. These improvements resulted in location of transition state geometry for 93% of cases and generation of kinetic parameters for 74% of cases. Newly calculated parameters agree well with benchmark calculations and perform well when used to replace estimated parameters in a detailed kinetic model of methanol combustion.


Introduction
Detailed kinetic models allow researchers to understand the chemistry of complex phenomena in systems such as combustion and hetrogeneous catalysis, thus enabling them to make informed experimental design choices, and to design and optimize processes and devices. Microkinetic models often contain hundreds of intermediates and thousands of reactions, for which thermodynamic and kinetic parameters need to be specified [1,2]. These parameters are ideally determined experimentally or calculated theoretically with high accuracy, but most are estimated [3]. These estimations allow parameters to be determined quickly, but usually with less fidelity [4]. Thermochemistry estimates are often derived from Benson's group additivity, where groups of atoms with known thermochemistries are summed [5]. These estimates are reasonable for most situations, but have been difficult to extend to some cases such as polycyclic species [6], motivating automated quantum mechanical or semi-emperical calculations [6,7].
When estimating the kinetics of a reaction, the Evans-Polanyi relationship can be used to estimate kinetics based on the enthalpy change of a specific reaction [8], if the rates of similar reactions are sufficiently well known. Alternatively, group contribution methods can provide kinetic estimations in a similar fashion to Benson's additivity methods [9][10][11]. Both Evans-Polanyi relationships and group additivity methods are fast and easily automated, which is especially useful in the generation of microkinetic mechanisms. Unfortunately, estimations fall short when exploring novel systems when the estimation rules are poorly known due to a lack of training data. In these cases, rule-based methods use less appropriate rules, and group-based methods utilize less specific group values, leading to errors in the rate as large as several orders of magnitude.
With growing computational power, calculating accurate kinetic parameters through transition state theory (TST) is no longer infeasible. However, TST calculations require a trained guess of transition state geometries, and often manual entry to arrive at trustworthy parameters. Given the number of reactions that are present in a detailed combustion model, these quantum calculations need to be automatized.
Automatizing TST calculations has been the focus of many [12]. This paper focuses on recent improvements to the AutoTST framework first developed by Bhoorasingh and co-workers [13,14].
AutoTST is an automated algorithm to locate reactant, product, and transition state (TS) geometries using quantum chemical calculations, to arrive at reaction rate expressions. AutoTST was originally built as a module within the Reaction Mechanism Generator (RMG) software [15] and could determine modified Arrhenius parameters from RMG reaction objects matched to one of three specific reaction families: unimolecular hydrogen migration (1 reactant to 1 product), radical addition to a multiple bond (2 reactants to 1 product), and bimolecular hydrogen abstraction (2 reactants to 2 products). Bhoorasingh noted that this was a sizable step in automated kinetic calculations, but the workflow needed some improvements. This work addresses these improvements to increase fidelity and speed of calculations, such as including a detailed conformer search, 1-D hindered rotor approximations, graph based symmetry number calculations, and parallelization of calculations using the Simple Linux Utility for Resource Management (SLURM) job scheduler. We observe these changes by recalculating reaction rates present in the Lawerence Livermore National Lab's (LLNL) butanol model by Sarathy and co-workers [16], comparing to benchmark calculations, and assessing the impact on iginition delay time predictions of the whole model.

Algorithm
The first generation of AutoTST is described by Bhoorasingh and co-workers [13,14]. This paper underscores important aspects of the original framework and highlight improvements. An overall updated workflow is described in Figure 1a with five broad steps.
First, initial geometry estimates are created (section 2.3) for both the reactants (using the distance-geometry and force-field methods in RDKit) and the TS (using the original AutoTST algorithm to modify the bounds matrix [13]). Then an ensemble of conformers are generated for each structure (section 2.4). The conformers are then processed in parallel to optimize the geometries (section 2.5) and validate them (section 2.6) to ensure the desired molecule or TS has been found.
Finally, temperature-dependent kinetics are calculated (section 2.7), including a correction to estimate the effect of hindered internal rotors. Figures 1b, 1c, and 1d are sub-workflows for the species, TS, and kinetics, that are performed during the overall workflow.

Inputs
AutoTST originally required users to provide reactions matched to one of three supported reaction families, or templates, present in RMG [15], and the electronic structure calculation settings resulted in long and complex input scripts. The user interface has now been updated, making inputs more straightforward. Users no longer have to provide matched, templated RMG reactions -the workflow will now automatically match the reaction of interest to one of the supported reaction families and identify the reacting atoms. Current supported reaction families are hydrogen abstraction, radical addition to multiple bonds, and intramolecular hydrogen migration ( Figure 2).
A user will supply a reaction of interest as an instance of an RMG Reaction class, or in a simple text string formed from the reactants and products in SMILES format (e.g.

CCO[O]+CC_CCOO+C[C]
for the reaction C 2 H 5 OO · + C 2 H 6 −− −− C 2 H 5 OOH + CH 3 CH 2 · ). If unable to match the input reaction to a supported reaction family, AutoTST will return an error. Reactions matched to   reaction templates are then used to create three dimensional TS geometries, described in subsequent sections.
Users are also required to specify electronic structure calculation settings, such as functional and basis set, and SLURM settings (e.g. username, partitions, accounts, excluded nodes) to enable AutoTST to perform calculations in parallel on computer clusters.

Representation of Species
The newest version of AutoTST uses the Atomic Simulation Environment (ASE) [17], RD-Kit [18], and RMG software packages to generate Species and Conformer objects that can be easily manipulated. A Species object requires a user to provide one or more SMILES strings. If multiple SMILES strings are provided, AutoTST will check that they are resonance isomers of each other; if only one is provided then additional resonance structures will be generated using RMG.
Once the list has been vetted, the Species object can generate Conformer objects for each resonance structure. Species objects act as a hierarchical class to organize one or multiple Conformer objects that represent the 3D geometry of a species.
To generate the 3D geometry of the conformer, AutoTST will first use the SMILES string to create an RMG molecule that describes the atom connectivity as a molecular graph. AutoTST then creates a 3D structure using the embed feature in RDKit. generate_conformers may be called to perform a systematic conformational analysis to create an ensemble of low energy conformer objects using an ASE electronic structure calculator. The conformer generation workflow is described in detail in the following sections.
Just as Species and Conformer objects are necessary for the representation of species and their conformers, AutoTST Reaction and TS objects are necessary to organize TSs and their conformers.
Reaction objects represent the chemical reaction, while the TS objects act as individual conformers for the reaction's transition state. Discussion of how TS geometries are constructed is described in the following section.

Creating Initial 3D Geometries
Initial species and TS geometries are generated as they were in Bhoorsingh and co-workers' original workflow [13,14]. For stable species, the embed feature in RDKit [18] is used to generate reasonable geometry guesses; TS geometries require additional treatment.
For TSs, reacting atoms are identified using the matched reaction family. In the case of hydrogen abstraction, these reacting atoms represent the abstracted hydrogen, the atom bound to the abstracting hydrogen, and the radical atom abstracting the hydrogen ( 2 H, 1 R, and 3 R· in Figure 2, respectively). The TS complex is then passed to a hierarchical decision tree that provides guesses of key distances between reacting atoms based on the reaction family and functional groups near the reaction center. The decision tree is descended to find the node where the functional group and its proximity to the reaction center closely matches the TS of interest. This node then provides the distances between the reacting atoms and serves as our "key distances".
A 3D geometry of the TS complex is created using the geometry embedding feature in RDKit [18] to generate a distance matrix. The distance matrix is a square matrix that describes the maximum and minimum allowable distances between pairs of atoms. The entries in the distance matrix are edited such that distances between reacting atoms are specified using the "key distances" determined by the decision tree. A new constrained TS geometry is generated by RDKit using the edited distance matrix.

Conformer Analysis
The original version of AutoTST would generate Species and TS geometries at random using RDKit [18] and would not ensure that AutoTST found the lowest energy, or most probable, conformation for each geometry [14]. AutoTST now performs a systematic conformer search on species and TSs by considering all dihedral angles not in a ring or containing a terminal methyl group, all invertible double bonds, and all chiral centers. Possible geometries are generated by creating all combinations of dihedral angles (range 0°to 300°with a 60°spacing), invertible double bonds (E vs Z configuration), and invertible chiral centers (R vs S configuration).
These initial geometries are optimized using Hotbit [19], a density-functional tight-binding calculator, and the BFGS optimizer provided in ASE [17] with a maximum of 1000 optimization steps.
Through this workflow, all optimizations are to energy minima where Species optimizations are completely unconstrained but TS geometries undergo a constrained optimization with the distances between reacting atoms fixed. Optimized conformers are compared against the initial geometry to ensure that isomorphism is maintained. Isomorphic conformers are then compared to identify unique conformers by asserting the average root mean square deviation between all other conformers is greater than 0.5 Å. All unique conformers within a specified energy cutoff of the lowest energy conformer are further optimized using DFT. For this study, we used 10 kcal/mol as our energy cutoff.

Optimize Geometries
Both versions of AutoTST relax species to minima but handle transition state geometry optimizations differently [14]. After determining the initial TS geometry from a group-estimation tree, the original workflow would perform three consecutive geometry optimizations on the TS. The first optimization froze distances between reacting atoms, while relaxing all others to an energy minimum. The second optimization froze distances between non-reacting atoms and relaxed all others to a saddle point. Lastly, the entire geometry is relaxed to a saddle point. To reduce computational costs, the new workflow skips the reaction center optimization.

Validate Geometries
The original AutoTST [14] would validate TS geometries by performing intrinsic reaction coordinate (IRC) calculations [20] and compare the output geometries against the input reactants and products. If the geometries matched, then the TS is validated.
These calculations were the bottleneck of the previous workflow and needed to be improved.
Inspired by the procedure of Van de Vijver in Genesys [21,22], rather than performing an IRC calculation on each TS identified, a vibrational vetting step is performed and followed by an IRC calculation if vetting is inconclusive. This vibrational vetting step entails reading in the log file of the saddle point and ensuring that there is: 1) only one negative frequency, 2) the distance between reacting atoms is less than two Å, and 3) the change in bond length by vibrational translation for reacting bonds is an order of magnitude greater than the change in bond lengths of non reacting bonds. If this vetting step is passed, no IRC calculations are performed because it is assumed the saddle point corresponds to the reaction of interest. If this vetting step is not passed, an IRC calculation is performed and used to validate the saddle point.

Kinetics Estimation
AutoTST uses the software package Arkane to estimate kinetic parameters [14]. Arkane is bundled with RMG [23], and can be used for pressure-dependent rate calculations [24] but here is used for canonical TST calculations Following traditional TST, it is assumed that the reactants and the TS are in a quasi-equilibrium state [25] in a vacuum and the rate limiting step is the transition from the TS to the products. From these assumptions, a modified form of the Eyring equation (equation 1) is used to relate the thermodynamic properties of the TS and reactants to the elementary rate of reaction [26]: where k B is the Boltzmann constant, T is the temperature, h is Plank's constant, R is the gas law constant, κ is the correction faction for quantum tunneling, and ∆G ‡ is the change in Gibb's energy between the TS and the reactants. Factors such as the change in Gibb's free energy and tunneling correction are conformer dependent and, as such, will be different in the new AutoTST workflow.

Symmetry Number Calculations
Symmetry numbers are a measure of the number of indistinguishable orientations a molecule or TS geometry can have [27]. Symmetry numbers were previously calculated using the SYMMETRY package which uses a 3D based approach [28]. However, symmetry numbers were often under estimated because the 3D comparison would often break symmetry for minor deviations. To combat this, we removed the SYMMETRY package from our workflow and now use the point group calculator included in Arkane. This allows Arkane to automatically calculate symmetry numbers when calculating kinetic and thermodynamic parameters, based on an analysis of the molecular graph or connectivity.

Hindered Rotor Correction
AutoTST previously used the rigid-rotor harmonic oscillator (RRHO) approximation when calculating kinetics, but this approximation can be inaccurate for geometries with internal rotors [14,24]. For most cases, the one-dimensional hindered rotor (1DHR) approximation is a more accurate representation of the internal rotation of stable species and TSs. We have added a modified 1DHR workflow using the approximation provided by Cohen [29] to account for internal rotation, without performing 1DHR scans as these are difficult to automate for saddle points. This correction is performed on kinetics after they have been calculated using Arkane [24].
First, we identify all rotatable torsions that are not present in a cycle, and for each rotor, we count the minimum number of substituents on either end of the rotor. This number is used to estimate the barrier to internal rotation using estimates from Benson [30] ( Table 1).
3.5 Table 1: Estimated barrier heights to internal rotation based on number of substituents.
The reduced internal moment of inertia and the internal symmetry of the torsion are calculated using methods available in RMG and these, with the barrier height, are used to calculate the approximate vibrational frequency using Eq. (2) [29] and the free rotor partition function using Eq. (3): (2) In Eq. (2) and Eq. (3), σ int is the internal symmetry of the torsion, I is the reduced internal moment of inertia, V is the barrier height, R is the gas law constant, k B is the Boltzmann constant, T is the temperature, and h is Planck's constant.
Values for vibrational frequencies, barrier heights, and free rotor partition functions of each rotor are used with Table 2 and Table 3 [29] to interpolate the vibrational and the hindered rotor contribution to the rate constant, k vib and k h.i.r. respectively.    Table 3: Best three-parameter fits for hindered internal rotor contribution to k(T ). From [29] For each rotor of each reactant and TS, the ratio of these numbers is calculated over temperatures ranging from 298-2500 K and is used to correct the rate constant calculated by the RRHO approximation: where k RRHO is the temperature dependent rate constant calculated using the RRHO approximation, k h.i.r,i is the 1DHR contribution to the rate constant for rotor i, k vib,i is the RRHO contribution to the rate constant for rotor i. For reactants, the product of these ratios over N reactants number of reactants is taken. Modified rates are fit to a three-parameter Arrhenius expression and returned to the user.

Methods
To test the efficacy of recent changes in AutoTST, Sarathy and coworker's model for the combustion of butanol was revisited [14,16]. To assess our changes, we wanted to study three distinct categories: 1. Success rate: how many reactions were we able to obtain kinetics for?
2. Micro-effects: how do AutoTST calculated rates compare to a set of benchmarks?
3. Macro-effects: how do AutoTST calculated rates impact an observable like ignition delay?
We attempted calculations on reactions present in the previous model used for the original AutoTST study, compared these calculations to benchmarks, and used these calculations to observe the change in ignition delay. These steps are described in detail in subsequent sections.

Success Rate
To assess the improvements made in AutoTST, we used the same set of 1117 reactions used in the original AutoTST paper plus three new hydrogen abstraction reactions that we were able to find using an updated version of RMG. These reactions come from Sarathy and coworker's model for the combustion of butanol [16] and are summarized in Table 4. Calculations were repeated on these reactions and used to observe improvements in our success rate in finding kinetics.

Model Chemistry
Reactions calculated from the original AutoTST paper [14] were performed using the M06-2X [31] functional with the MG3S [32] basis set. For this study, calculations were performed using the using the M06-2X function with the cc-pVTZ [33] basis set. Different model chemistry was used because Arkane no longer supports the MG3S basis set.

Benchmark calculations
As part of the validation methods from the original AutoTST paper, Bhoorsingh and co-workers performed a series of benchmark calculations on six reactions shown in Table 5 [14]. These reactions represent two reactions from each of the three supported reactions families where the AutoTST rate expression disagreed with RMG predicted expressions by more than a factor of 100 at 1000 K and 1 bar.

Label Family Reaction R1
Hydrogen Abstraction Hydrogen Abstraction These benchmark calculations were performed using the original TS calculated from the original AutoTST workflow, but with 1-D hindered rotor calculations performed on each rotatable dihedral.
If a lower energy conformer was identified during a scan, the lower energy conformer was optimized and 1-D hindered rotor scans were restarted on the new conformer.
Single point energies were calculated using ORCA [34,35] at the CCSD(T)-F12/RI method with the cc-VTZ-F12 [36] and the cc-VTZ-F12-CABS [37] basis sets. In addition, point groups for symmetry numbers were determined by hand. These additional treatments were used to calculate rate expressions with Arkane and are referred to as the "Benchmark" calculations. In this work, these benchmark calculations were compared to the original AutoTST, LLNL, and new AutoTST kinetic parameters and discrepancies were measured.

Generation of alternative models
To observe the impact of AutoTST calculations, we generated alternate models using our newly calculated AutoTST kinetics. Kinetics that were successfully calculated by the improved workflow were swapped into the LLNL butanol model [16] one at a time to generate alternate models. E.g.
if we calculated kinetics for reaction X, we generated an alternative model by swapping in kinetics calculated by the original workflow and the new workflow. These models were tested against four sets of ignition delay data using PyTeCK [38] (described in the following section) to quantify error of our alternative models against experiments.

PyTeCK
PyTeCK, Python tool for Testing Chemical Kinetics [38], was used to quantify the error in theoretical models against experimental data. PyTeCK works by reading in experimental data in the human-and machine-readable ChemKED, Chemical Kinetic Experimental Data, format [39].
PyTeCK will read the experimental conditions from a ChemKED file, perform the simulation using Cantera [40] at the experimental condition, and return the error between the simulated value and the experimental value as described in equations 5 and 6. PyTeCK is currently limited to ignition delay experiments.
In equations 5 and 6, E j is the average error for the jth data set with n data points, τ exp i is the experimental ignition delay measured at conditions i, τ sim i is the simulated ignition delay at conditions i, σ exp i is the uncertainty of the the experimental measurement for the ith data point, E tot is the total error over m data sets.
PyTeCK was used to compare alternative models against four sets of experimental summarized in Table 6.

Workflow Efficacy
One goal of this work was to improve the success rate of AutoTST, which is measured here in two ways: (1) observe the percentage of reactions where AutoTST found a TS and (2) observe the percentage of reactions where AutoTST arrived at a rate expression. These results are summarized in Table 7 and   A different trend is noticeable when observing the percentage of reactions where AutoTST found a rate expression. For our re-run, calculations were attempted on 1120 reactions in the LLNL butanol model [16] and rate expressions were found for 833 reactions or 73% -a small increase compared to the original workflow. In the original workflow, Bhoorasingh and co-workers found the success rate was independent of reaction family [14], but here, there is an inverse correlation between the success rate and the number of reactants and products. Hydrogen abstraction reactions (two reactants and two products) had a success rate of 70.1%, radical addition reactions (two reactants and one product) had a success rate of 83.2%, and hydrogen migration reactions (one reactant and one product) had a success rate of 92.3%.
There were 208 failures when calculating kinetics from an apparently successful TS optimization.
197 of these failures were barrier height errors (186 from hydrogen abstraction and 11 from radical addition reactions). These failures occur when either the forward or reverse barriers in the Eckart model are negative, indicating that the energies of the reactants or products is greater than the TS. This can be because the TS complexes connects to two van der Waals (vdW) wells rather than the bimolecular entrance or exit channels. These wells allow for intermolecular interactions to occur between reacting species that lower the potential energy of the complex, leading to a saddle point that is submerged below the reactants or products. These errors were more common in abstraction reactions with oxygen-containing reactants (e.g abstraction by O --O, • OH, • OOH) which have strong vdW interactions. These errors might be addressed by treatment of vdW wells with a 3TS model solved with a Master Equation (ME), as is done in EStokTP [46]. These barrier height errors account for almost all of the failures in the hydrogen abstraction reaction family.
The remaining 11 errors are attributed log parsing errors that occur when AutoTST interfaces with Arkane. These errors can be remedied through development of Arkane, RMG, and AutoTST in tandem.

Benchmark Calculations
We revisited the six benchmark calculations that were performed in the original AutoTST study [14]. This set of reactions consists of two representatives from each reaction family supported by AutoTST. For each of these reactions, the rate coefficients estimated by RMG and from the original AutoTST disagreed by an order of at least 1000 when calculated at 1000 K. Arrhenius plots were generated using rate expressions from this work, the original AutoTST study, the original butanol model, and the benchmarks from the original AutoTST study in Figure 3. From Figure 3 we see good agreement between newly calculated rate expressions and the bench-marks. To quantify the agreement, a discrepancy measurement was defined as: where d i (T ) is the discrepancy between the benchmark rate expression, k 0 , and the rate expression for the ith source, k i , at temperature T . This discrepancy represents an order of magnitude difference between a rate expression and the benchmark. The average difference between the two rate expressions is calculated over temperatures between 300 K and 2000 K with 100 K spacing. Table 9 shows the average and standard deviations in discrepancies between kinetics from the original workflow, the updated workflow, and the LLNL butanol model. A high standard deviation indicates a difference in apparent activation energy.  Table 9: A table containing the discrepancies between kinetic parameters and benchmark calculations from [14].
In five out of six benchmark reactions, the updated workflow had the lowest average discrepancy overall. Reaction 2 was the only exception where rate coefficients from the original workflow performed the best. In all cases, the updated workflow performed better than the LLNL butanol model. When taking into account the variability of the discrepancy measurement, it is difficult to conclude that any one set of kinetics performed the best. However, it is clear that the modifications have resulted in more accurate kinetics as a whole.

Comparison against experimental data
We generated 833 alternate models for the combustion of butanol by swapping in our newly These tables also summarize the absolute difference in log 10 k at 1000 K and 100 bar. In addition, reactions with the greatest disagreement in log 10 k are reported in Table 12.   In all reactions present in Table 10, rate rules or approximations were replaced by calculated rates. It is reasonable to deduce that the reduction in error was most likely attributed to this change in source. The "Rate Rules" entry in LLNL Source column refers to a 1998 rate rule developed for low temperature isomerization of RO 2 species in an n-heptane microkinetic model [47]. This rule creates a two-parameter Arrhenius expression by using the change in enthalpy for the reaction, a  pre-exponential factors. The "Approximation" entry indicates that parameters come from educated approximations by the creators of the model.
It is also understandable that AutoTST calculated parameters increased error in some cases.
LLNL kinetics from four out of the five reactions present in Table 11 come "Computational Studies" [48][49][50]. The LLNL parameters from "Computational Studies" were TST calculations that were created and analyzed by hand. It is prudent to say that kinetic parameters from a thorough computational study would best any automated rate calculation; therefore, the increase in error upon using the AutoTST calculations is reasonable. With regards to the "Database" kinetics, these ultimately led to a 1972 book that is out of print [54].
Finally, when investigating the sources where kinetics disagree the most, another trend is seen.
Two out of five of the LLNL kinetics come from approximations, with another two coming from "Analogies" (i.e. modification of existing parameters using a trained guess). Between approximated and calculated parameters, large differences in kinetics are to be expected; although the analogies between kinetics of different sources can be reasonable estimates, they are still estimates so it is understandable why the kinetics disagree so much. Interestingly, we were unable to trace back the source to LLNL reaction 1489 completely. After a detailed literature search, this number originated from a kinetic model for the combustion of hexane [52], but found its way to the butanol model through a reference of a reference. Even though these kinetics disagree by 10 or more orders of magnitude, the changes did not change the PyTeCK error by much, indicating that the ignition delay is not sensitive to these reactions.
These results highlight that kinetics calculated from AutoTST perform as expected when utilized in a detailed kinetic model -they perform poorly in comparison to parameters from careful computational studies, but improve upon kinetics that come from estimates or approximations.
They also improve upon calculations from the first version of AutoTST.

Conclusions
Through this work, improvements have been made to AutoTST and include: rewriting of the code to make user inputs simpler (as well as developing, debugging, and testing), inclusion of a systematic conformer search, modification of how symmetry numbers are calculated, and addition a 1-D hindered rotor treatment. These changes were measured by attempting automated calculations on 1120 reactions present in the LLNL model for the combustion of butanol [16], the model investigated in the original AutoTST study [14].
Of the 1120 reactions, 858 were bimolecular hydrogen abstraction reactions, 78 were unimolecular hydrogen migration reactions, 184 were radical addition reactions. The updated version of AutoTST was able to find TS geometries for 92.9%, or 1041 reactions: 797 (93%) for hydrogen abstraction, 73 (94%) for hydrogen migration, and 171 (93%) for radical addition. The updated version of AutoTST was able to arrive at kinetic parameters for 608 (70.9%) hydrogen abstraction reactions, 72 (92.3%) hydrogen migration reactions, and 153 (83.2%) radical addition reactions. This dependence of success rate on reaction family is most likely tied to the relation of the energy difference between the reactants, products, and TS geometry, with abstraction reactions more often passing through van der Waals wells leading to a saddle point submerged below the reactants. The overall success rate was an improvement in comparison to the original version of AutoTST.
Six benchmark reactions identified in the original AutoTST study were compared against parameters from the LLNL model, the original AutoTST calculated kinetics, and the updated kinetics.
For five out of six of these reactions, the updated version of AutoTST had the lowest discrepancy (i.e. agreed the most with the benchmarks) and for all six reactions, the updated version of Au-toTST had a lower discrepancy than LLNL parameters. This shows that these modifications to AutoTST resulted in more accurate kinetics that are closer to "Chemical Truth." We also observed the effect of our calculated parameters on ignition delay. By utilizing PyTeCK, we saw the impact of our calculations when implemented in a detailed kinetic model. In 94% of cases our modified kinetic parameters had a negligible impact on ignition delay when applied individually.
AutoTST was able to decrease the error for reactions that were sourced from either approximations or rate rules, but increased the error for reactions that came from computational studies, which is to be expected. AutoTST is designed to efficiently provide on-the-fly kinetic parameters that perform better than estimates, but will most likely not outperform thorough computational studies.
When comparing LLNL and updated AutoTST kinetics there are some kinetic parameters that disagreed by over 10 orders of magnitude. In these cases the LLNL parameters originated from approximations or analogies so it is likely that the parameters from AutoTST are more accurate.
However, they have little impact on ignition delays for this system.
Although many improvements have been made to AutoTST, there are still more to be made.
The systematic conformer search is adequate for small species and TSs but is too computationally expensive for larger complexes, so stochastic methods should be investigated. In addition, the hindered rotor treatment in AutoTST could be improved. Proper rotor scans should be performed to accurately account for hindered rotor effects in both species and saddle points. High fidelity singlepoint energy calculations could increase accuracy. Finally, treatment of vdW wells and multiple TSs would improve the calculation of abstraction reactions.
In addition to the developments of AutoTST, there are a number of future research directions that can make use of AutoTST or automated rate calculators. AutoTST can be used to help researchers identify and resolve errors in detailed kinetic models by (1) finding discrepancies in kinetic parameters and (2) utilizing updated parameters in a detailed kinetic model to bring it closer to the Chemical Truth. Finally, using AutoTST in tandem with an automated reaction mechanism generator like RMG to create kinetic models, where a majority of rates are currently estimated, would be fruitful.