A General User-Friendly Tool for Kinetic Calculations of Multi-Step Reactions within the Virtual Multifrequency Spectrometer Project

We discuss the implementation of a computer program for accurate calculation of the kinetics of chemical reactions integrated in the user-friendly, multi-purpose Virtual Multifrequency Spectrometer tool. The program is based on the ab initio modeling of the involved molecular species, the adoption of transition-state theory for each elementary step of the reaction, and the use of a master-equation approach accounting for the complete reaction scheme. Some features of the software are illustrated through examples including the interconversion reaction of hydroxyacetone and 2-hydroxypropanal and the production of HCN and HNC from vinyl cyanide.


Introduction
Calculation of chemical kinetics in the gas phase by accurate theoretical models is extremely important in research areas like atmospheric chemistry, combustion chemistry, and astrochemistry. As a matter of fact, the accurate prediction of reaction rates and the evolution of the involved species in a given set of physical conditions is a key feature for understanding the presence of a molecular or ionic species in that environment. Sometimes, the reactions involved are too fast to be tracked by laboratory experiment or the associated physical conditions are simply not reproducible, hence the understanding of those reactions relies on accurate theoretical treatments capable of predicting the evolution of the species involved in a reaction network.
A rigorous treatment of the time evolution of a chemical reaction should be based on quantum-dynamics calculations modeling the motion of the involved nuclei on ab initio potential-energy surfaces [1][2][3]. However, exact quantum-dynamics methods are only applicable to very small systems made up of three or four atoms (see, for instance, Refs. [4][5][6]). For the remaining systems-the vast majority-one can either opt for approximate methods, such as the Multi-Configuration Time-Dependent Hartree (MCTDH) [7] or the Ring Polymer Molecular Dynamics [8] (which can however extend this limit only marginally), or for a (quasi-)classical treatment of the nuclear motion [9][10][11].
On the other hand, a less expensive yet reliable route to chemical kinetics is the adoption of statistical models, such as the popular transition-state theory (TST) in one of its variants, which successfully exploits information on the energetics of a limited set of important points of the potential energy surface to predict the kinetics of chemical reactions. The usual procedure in this framework involves the calculation of transition states and intermediates of a given reaction and a description of the motions at molecular level of these species. Then, classical or semiclassical transition state theory (TST) is applied to calculate the reaction rates of each of the elementary steps making up the whole reaction (the Rice-Ramsperger-Kassel-Marcus (RRKM) [12][13][14] theory, shortly summarized in the following, is usually adopted for unimolecular reactions in the gas-phase). Finally, the time evolution of the relative abundances of each of the reactant, intermediate, and product species is calculated using methods based on either master-equation or stochastic approaches.
In this paper, we discuss the implementation of the computer program StarRate for kinetics calculations of multi-step chemical reactions, and its integration in the graphical interface of the user-friendly, multipurpose framework Virtual Multifrequency Spectrometer (VMS) [15]. The Virtual Multifrequency Spectrometer (VMS) is a tool that aims at integrating a wide range of computational and experimental spectroscopic techniques with the final goal of disclosing the static and dynamic physical-chemical properties of molecular systems, and is composed of two parts: VMS-Comp, which provides access to the latest developments in the field of computational spectroscopy [16,17], and VMS-Draw, which provides a powerful graphical user interface (GUI) for an intuitive interpretation of theoretical outcomes and a direct comparison to experiment [18]. We discuss the integration of StarRate within the VMS tool and illustrate some features of the developed software through two important reactions: the single-step interconversion of hydroxyacetone and 2-hydroxypropanal, and the more challenging multi-step dissociation of vinyl cyanide. It is worth mentioning here that the reported calculations were performed for the purpose of illustrating the developed computational software, and that providing new accurate results on the above reactions for comparison with experiment is beyond the scope of this work.
The article is organized as follows. Section 2 is devoted to computational details of the developed software. In Sections 3 and 4, we address the kinetics of the above mentioned reactions. In Section 5, conclusions are drawn and perspectives are outlined.

Computational Details: StarRate and the VMS Tool
StarRate is an object-based, modern Fortran program for modeling the kinetics of multistep reactions. At its current stage of development, StarRate targets multi-step unimolecular reactions (which can however dissociate, irreversibly, to multiple products). (The implementation of procedures for the treatment of bimolecular entrance channels is currently in progress). From a technical point of view, the program is written in the so-called 'F language' [19,20], a carefully crafted subset of Fortran 95, and is conceived in an object-based programming paradigm. As described in deeper detail in Ref. [21], StarRate is structured in three main modules, namely molecules, steps and reactions, which reflect the entities associated with a multi-step chemical reaction. All of these modules contain a defined data-type and some related procedures to access and operate on it. The main program, StarRate, controls the sequences of the calling of the procedures contained in each of the three main modules.
Another important module of StarRate is in_out, which handles the input and output operations of the program. Input data are accessed by StarRate through an XML interface based on the same versatile hierarchical data structure that is adopted by VMS. (The current version of VMS reads data in the JSON format, so that a straightforward conversion from XML to JSON is in order. This can be easily done, for instance, using xml2json (https://github.com/hay/xml2json), or online converters such as https://www.convertjson.com/xml-to-json.htm). At the beginning, the user has to prepare a very simple input file encoding a reaction scheme (see the starrate.inp box in Figure 1) and gather all the files, one for each molecular species, containing data deriving from electronic-structure calculations. These can be either in an internal standard format (similar to that adopted in the EStokTP [22] package) or directly output files of quantum-chemistry packages, as exemplified in Figure 1 for the case of the Gaussian package. Currently StarRate supports output files from this quantum-chemistry package (.log extension), though support for other popular electronic-structure programs is presently being pursued (see also Refs. [23,24] on the issue of interoperability and common data formats in quantum chemistry). Then, a Python script is run which extracts data from the output files generated by the quantum-chemistry calculations and, driven by the reaction scheme, collects the information in the proper sections of the XML file. The structure of the XML interface is schematized in Figure 2. The whole XML document develops under a root element named escdata. The escdata has one child element for each molecule named section_run. Each of these elements contains three nodes: program_info, section_system, and system_single_configuration_calculation. All the information regarding a molecule is handled by these three sibling nodes. The program_info node contains two subnodes which keep track of quantum chemical software name and .log file location. The section_system node contains basic information which does not require quantum chemical information (viz., molecular charge, spin multiplicity, atom label, atomic numbers, rotational constants). The last sibling, system_single_configuration_calculation, contains information which requires quantum chemical calculations (viz., vibrational constants, SCF energy, density of state data). Finally, the last section_run collects information on physical conditions and on the reaction under study. For illustrative purposes, the actual .xml file for the reaction studied in Section 3 is given as Supplementary Material.
Once the XML has been generated, the StarRate program comes into play. The module in_out reads the XML file (a well-built external Fortran library, FoX_dom [25], is used for XML parsing) and saves the data for each molecule and step as structured arrays of the molecules and steps modules, respectively. Some information, such as vibrational frequencies, rotational constants, and electronic energy, is collected from the electronic-structure calculations; some other information such as densities of states (see later on) and single-step microcanonical rate coefficients are either also read as input data or computed internally to StarRate. Lastly, the reactions module solves the kinetics for the overall reactions using a chemical master equation method. At the end of the calculations, VMS is used to access, visualize and analyze the data produced thanks to the shared XML interface (see Figure 1).

Elementary Steps: The Interconversion Reaction of Hydroxyacetone and 2-Hydroxypropanal
The interconversion reaction between hydroxyacetone and 2-hydroxypropanal is an important reaction in the context of atmospheric chemistry because the hydroxyacetone represent the simplest form of photochemically oxidised volatile organic compounds [26]. In a recent study, Sun et al. [27] have considered the interconversion mechanisms on several hydroxycarbonyl compounds, and much attention has been focused on the interconversion reaction between hydroxyacetone and 2-hydroxypropanal. This isomerization reaction can occur through three different mechanisms, 2 high-barrier multistep processes and, as shown in Figure 3, a direct mechanism via double hydride shift involving a low-barrier concerted transition state. In their work, Sun et al. also supposed that hydroxycarbonyl compounds can adsorb solar radiation, as carbonyl compounds, from 320 to 220 nm and then undergo an internal conversion to the vibrationally excited ground state with an energy more than sufficient to overcome the isomerization barrier, and computed RRKM microcanonical rate coefficients in order to understand how much the isomerization reaction is favored with respect to collisional deactivation and fragmentation processes at a given excitation energy.
Within the RRKM theory [12][13][14], the microcanonical rate coefficient for the reaction of Figure 3 is given by the equation [28] k where In Equations (1) and (2), h is Planck's constants, N ‡ (E) is the sum of states of the transition state (TS) (computed by excluding the normal mode with imaginary frequency under the assumption that the motion along the reaction coordinates is separable from that of the other modes), and ρ(E) and ρ ‡ (E) are the density of states (DOS, i.e., the number of rovibrational states per energy interval) of the reactant molecule and transition state, respectively. As apparent, a central quantity in this framework is the molecular rovibrational density of states of the involved molecular species. This can be easily worked out by convoluting its rotational and vibrational counterparts [29]. In the present version of the program, a classical expression is used for the rotational DOS (see [21] for details), while the vibrational DOS is evaluated at uncoupled anharmonic level by adoption of the Stein-Rabinovitch [30] extension of the Beyer-Swinehart algorithm [31]. An improved version of Equation (1) accounts for the tunneling correction by using a modified version of the sum of states N ‡ (E) of the TS. A common and efficient way of including tunneling is by means of the asymmetric Eckart barrier [32]. Within this model, the sum of states of the transition state is redefined by where N ‡ tunn (E) is a tunneling-corrected version of the sum of state of the TS and V 0 is the classical energy barrier for the forward reaction. The quantity P tunn (E ) is the tunneling coefficient at the energy E , and is given by the expression where, a, b, and c are parameters defined by: Here, V 1 is the classical energy barrier for the reverse reaction, and ν i is the magnitude of the imaginary frequency of the saddle point (in Equation (5), h = 1 if the energies are expressed, as in this work, in cm −1 ). For illustrative purposes, we computed the microcanonical rate coefficient for the direct and inverse reaction of Figure 3, both with and without tunneling correction. To this purpose, the three molecular species were modeled by density-functional theory with the B2PLYP-D3/jun-cc-pVTZ model chemistry. According to our calculations, the forward reaction is exoergic by 1719 cm −1 with a barrier of 16448 cm −1 (the barrier for the backward reaction is 14729 cm −1 ). The resulting microcanonical rate coefficients k(E) are plotted in Figure 4 (on a logarithmic scale) for the forward (blue line) and backward (red line) reaction as a function of the energy relative to the hydroxyacetone zero-point energy. In the same figure, the dashed curves are the tunneling-corrected ones. As apparent, the tunneling correction enhances the reaction rate both in the forward and backward direction, more visibly nearby the threshold region, thus lowering the actual value of the reaction threshold in both directions. The thermal rate coefficient can easily be computed from the microcanonical rate coefficients using the following equation: with Q(T) being the partition function of the reactant species. The computed thermal rate coefficients for the forward reaction (isomerization reaction of hydroxyacetone to 2-hydroxypropanal) in the temperature range 151-501 K are given as Arrhenius plot (log 10 k(T) versus 1/T) in Figure 5.
Results are reported both neglecting (blue triangles) and including (blue circles) tunneling. These data can be fitted to the popular Arrhenius equation: (with A being the pre-exponential factor, E a the activation energy, and R the gas constant) or to the more refined Arrhenius-Kooij formula [33] (also known as modified Arrhenius equation [34]) allowing for a temperature dependence of the pre-exponential factor: which essentially implies a linear variation of the activation energy with the temperature, E a /R = γ + βT. The Arrhenius best-fitting curve for both the tunneling-corrected and no-tunneling data are shown in Figure 5 as dashed black line and solid black line, respectively. The Arrhenius-Kooij best-fitting curve for the tunneling-corrected data is also reported as a red dashed line in the same figure, while the best-fitting parameters together with the associated coefficient of determination R 2 are given in Table 1.

Arrhenius-Kooij (Equation (8)), with Tunneling
As evident from Figure 5 and Table 1, the Arrhenius equation perfectly fits the thermal rate coefficients calculated by neglecting tunneling, yielding a R 2 = 1.0000 and an activation energy of 16428 cm −1 that compares well with the computed reaction barrier. On the contrary, the tunneling-corrected thermal rate coefficients show a deviation from linearity with decreasing temperatures. As a result, the Arrhenius expression yields a worse best-fitting curve (R 2 = 0.9973) and a lower activation energy of 15120 cm −1 , while a better fitting is obtained through the Arrhenius-Kooij equation (R 2 = 9.9992), which gives an activation energy of 14839 cm −1 at T = 250 K and of 13352 cm −1 at T = 150 K.

Multi-Step Reactions: The Dissociation of Vinyl Cyanide
The dissociation of vinyl cyanide (VC, C 3 H 3 N), is particularly interesting because it involves multiple reaction channels and different sets of products (HCN, HNC, HCCH, and :CCH 2 ) and hence it serves as a very good test case for master-equation based kinetic models. The potential-energy surface for this reaction has been investigated in a recent work by Homayoon et al. [35] through ab initio CCSD and CCSD(T) calculations with the 6-311+G(2d,2p) and 6-311++G(3df,3pd) basis sets. In the same work, a reaction scheme involving ten unimolecular steps, three of which reversible, was proposed. The ten reaction steps are summarized in Table 2, while the associated reaction diagram is given in Figure 6. Table 2. Reaction steps involved in the dissociation mechanism of vinyl cyanide considered in this work. The associated activation energies (relative to the zero-point energy of the reactant of each step) are also given.

Reaction
Step As shown in Figure 6 and Table 2, VC can directly dissociate to product sets :CCH 2 + HCN and HCCH + HCN through steps 1 and 2 (the only direct dissociation paths), or lead to formation of reaction intermediates Int1-III (the most stable one), Int1-IV, and Int1-V, further evolving to products. On the other hand, product HNC can only be formed via indirect dissociation paths involving the above-mentioned intermediates. A screenshot of VMS showing the structures of all the molecular species involved in this reaction is given in Figure 7.
Within a master-equation approach (see for instance [36]), to determine the time evolution of the relative abundance of the involved species, initially a matrix, K, is set up by opportunely combining the microcanonical rate coefficients at a specified energy. In particular, the diagonal elements K ii contain the loss rate of species i, while the off-diagonal elements K ij contain the rate of formation of species i from species j. The rate of change in the concentration of each species is given by the vector differential equation: where c is the vector of the concentrations of the species at time t. This is a linear differential equation and can be solved by diagonalization of K. In terms of the eigenvector matrix Z and eigenvalue vector Λ, the solution of Equation (9) reads: where c(0) is the concentration vector at t = 0. In this model, a fundamental hypothesis is that collisional relaxation occurs on time scales much shorter than those that characterize phenomenological kinetics [37]. It is worth mentioning here that a more general version of the master equation would involve diagonalizing a much larger matrix explicitly including collisional relaxation [36]. However, if the above-mentioned hypothesis holds, the resulting eigenvalues would appear in two separated sets: one made up by so-called internal energy relaxation eigenvalues (IEREs) and one made up by so-called chemically significant eigenvalues (CSEs). These last eigenvalues that relate to the phenomenological kinetics of interest in interstellar space and atmospheric studies would be identical to those obtained by solving Equation (9).   By using the methodology described, we computed the evolution of species with respect to time through the StarRate program using the structural parameters of the species given in Ref. [35], and computing the microcanonical rate coefficients through Equation (1). For the reader's convenience, we give the full form of the matrix K for this reaction (please note that expressions in square brackets, though spanning several rows, relate to single matrix elements and are shown as such to give a compact picture of the matrix): In our calculations, the initial concentration of VC was taken as 1.0 and the concentration of other species was set to 0.0. The relative abundances of the involved species as a function of time are plotted for two energies, namely E = 51764 cm −1 and E = 62000 cm −1 , in Figures 8 and 9, respectively. Int1III  Int1IV  Int1V  Int2IV  Int2V   TS1III  TS1II  TS1IV  TS1I  TS1VII  TS1VI   TS1V  TS2III  TS2IV  By inspection of the plots, a first remark is that there is a sudden spike of the concentration for Int1-III in a very small time range for both energies. This is because Step 3 involves a low activation energy (20076 cm −1 ) compared to Int1-IV (37494 cm −1 ) and Int1-V (36724 cm −1 ), and the intermediate Int1-III is relatively stable compared to the other two (the stability of Int1-III is also reflected by the long sigmoidal tail of the plot). At the considered energies, these last two species are virtually never present and as soon as formed evolve into products. The reaction paths involving these two intermediates (the only ones leading to formation of HNC) become increasingly important at E = 62000 cm −1 ; in fact, while the branching ratio HCN/HNC tends to a value of about 2.5 at E = 51764 cm −1 , a branching ratio of 1.9 is obtained at E = 62000 cm −1 . The predicted branching ratio at E = 51764 cm −1 nicely compares with the experimental estimate of 3.3 of Ref. [38] and the theoretical one of 1.9 of Ref. [35] (where, however, only vibrational densities and sum of states were taken into account in the calculation of the microcanonical rate coefficients), substantially improving over former theoretical calculations yielding a branching ratio of over 120 [39].

VC
An interesting feature offered by VMS is that of visualizing matrices in a 'heat-map' fashion through a color palette reflecting the actual value of the elements. A heat-map of K for the rate coefficient at E = 62000 cm −1 is given in Figure 10. Looking at the first column and recalling the meaning of the K ij elements, one can see that VC is directly converted to all the remaining species except HNC. The highest conversion rate (darkest gray) is towards Int1-III, while the rate of formation of the other two intermediates from VC is slower due to higher activation energies. Direct conversion to product sets containing HCN is also fast. On the contrary, as already mentioned, HNC is not formed directly from VC (blank square in position 7,1). Looking at the whole matrix, the darkest squares are those of the matrix elements connecting Int1-IV and Int1-V to products, due to the associated lowest activation energy. This is in line with the fact that, as also shown by Figure 8, these intermediates evolve rapidly to products. 14 Figure 10. Color map of the transition K matrix at energy E = 62000 cm −1 relative to the reactant zero-point energy as visualized in the VMS software. For each element of the matrix, the value of | log 10 K ij | in s −1 is plotted according to the gray scale on the right-hand-side of the plot area.

Conclusions
In this paper, the implementation of a computer program for chemical kinetics of multi-step reactions and its integration with the graphical interface of the Virtual Multifrequency Spectrometer has been discussed. The developed computational machinery is built around an input/output interface using a hierarchical data structure based on the XML language and shared in common with VMS. Details are given on the implementation of the calculation of microcanonical rate coefficients for the single steps of a unimolecular reaction, and on the modeling of the time evolution of the relative abundance of the involved species. The main features of the program have been illustrated through two example reactions, namely the atmospherically relevant interconversion between hydroxyacetone and 2-hydroxypropanal, and the production of HCN and HNC by dissociation of vinyl cyanide. Work is ongoing in our laboratory to account for bimolecular entrance channels, enhance the potentialities of the program, and integrate it in virtual-reality environments [40,41]. Funding: The research leading to these results has received funding from Scuola Normale Superiore through project "DIVE: Development of Immersive approaches for the analysis of chemical bonding through Virtual-reality Environments" (SNS18_B_RAMPINO) and program "Finanziamento a supporto della ricerca di base" (SNS_RB_RAMPINO).