First Principles Based Automated Kinetic Model Generation Using on-the-Fly Ab Initio Calculations

Introduction

 

Many kinetic models for combustion and pyrolysis processes have been published in the past decade, and a common challenge in most work is often being stressed, i.e. the shortage of accurate rate coefficients for a large number of elementary steps. Currently, many rate coefficients are calculated based on rough approximation methods, and their accuracy can be questioned when these approximations are used outside of the original application domain of the applied method. Furthermore, building these approximation methods requires a complete training and test set of reactions of which the rate coefficients are known from either experiments or theory. Manually obtaining all these rate coefficients is time consuming, and this work is aiming at automatically calculating a large number of rate coefficients using quantum chemistry, minimizing the manual work. The developed methods can furthermore be used after the model generation, i.e. when not all the kinetics are sufficiently accurate.

This work has been performed in the framework of Genesys¹, an automatic kinetic model generation tool developed at the Laboratory for Chemical Technology which employs both the rule-based and the rate-based termination criteria. The latter can in particular benefit from the automatic ab initio calculations in the search for more complete, more accurate and predictive kinetic models.

Methodology

 

Genesys offers a representation of chemical reactions, supported by the Chemistry Development Kit². Molecules, radicals and transition states are treated as mathematical graphs, in which the atoms are nodes and the bonds are edges of the graph. This is a very powerful molecular representation to convert reactants into products, to find sub-molecular patterns, to search in databases for appropriate data, etc. However, it is insufficient to perform ab initio calculations, where the coordinates of the atoms in the 3D space are required. The conversion of connectivity information to 3D coordinates is already well established for molecules and radicals, for which rule-based heuristic methods and numerical methods exist. “Distance geometry” is a method belonging to the latter group in which the coordinates are estimated by minimizing an error function describing all the interatomic distances compared to expected distances, which are tabulated. Additional optimization steps, among others using the Merck Molecular Force Field (MMFF)³, allow the generate a meaningful set of 3D coordinates, directly applicable for ab initio calculations. This has been implemented in Genesys.

To obtain accurate rate coefficients, quantum chemical calculations of transition states are also required. The 3D coordinates of atoms in transition states are more complex to deduce and cannot be summarized in the same simple rules used for molecules and radicals. Since Genesys has been programmed to offer the end-user flexibility in terms of which reaction families to use and how to constrain them, this philosophy is further employed for the generation of 3D coordinates for transition states. The end-user has to possibility to guess the interatomic distances of the atoms in the transition state that change in connectivity throughout the reaction, called the reactive atoms. This information, the so-called “template” is used in Genesys to obtain 3D coordinates of all the atoms via “distance geometry”, i.e. similar to stable molecules and radicals. The template can be built from a small set of manually calculated transition states or from the vast literature data on quantum chemical calculations. Refining the “distance geometry” output to a more meaningful set of coordinates is also done. For this, the MMFF is used in which equilibrium distances for the reactive atoms originate from the template, and the force constants are taken from analogous stable molecules and radicals.

Once initial 3D coordinates are available for molecules, radicals and transition states, quantum chemical calculations are performed using the Gaussian09 suite of programs⁴ in several steps. First, the geometry is energy-optimized using fast semi-empirical calculations such as PM3. Using the resulting geometry, the conformational space is scanned by rotating around single bonds in a systematic and exhaustive manner. Each of these conformers are optimized using DFT methods, and the geometries belonging to the lowest 5 kJ mol^-1 are used for high-level optimization, i.e. CBS-QB3 calculations. Finally, hindered rotor potentials are obtained to correct for the harmonic oscillator approximation.

For transition states, it is necessary to verify if the imaginary frequency indeed corresponds to the expected reaction. This is done in two steps. First, Genesys controls if the final geometry is close to the user-provided template. Second, the normal mode corresponding to the imaginary frequency is analyzed and it is evaluated whether the bonds that are broken or formed have a large contribution to this normal mode. For this, the normal mode in Cartesian coordinates is transformed into a normal mode in internal coordinates.

Via statistical thermodynamics, standard enthalpies of formation, standard entropies and heat capacities of all the molecules, radicals and transition states are calculated. The frequencies that resemble internal rotations are automatically searched for via an analysis of the normal mode. Their contributions to the vibrational partition function are automatically replaced by the hindered rotor partition functions. Also, the total rotational symmetry number of the species is automatically calculated and used to obtain correct entropy values. In the case where stereocenters are present, several approaches are possible. In Genesys, it is possible to explicitly track each stereoisomer as a separate species and to calculate separate thermodynamic data. However, when this is not necessary, they can also be lumped together, for which the entropy is then corrected by the number of optical isomers.

The thermodynamic data of reactants and the transition state of a reactions are subsequently used to obtain high-pressure limit rate coefficients. Tunneling corrections are automatically calculated using Eckart potentials, and Arrhenius parameters are determined by regression of the rate coefficients.

Results

 

The new methods have been illustrated by calculating rate coefficients for 8 reactions, each belonging to a different reaction family, i.e. inter- and intramolecular hydrogen abstractions, inter- and intramolecular radical additions, a substitution, a dehydration, a dehydrogenation and a retro-ene reaction. For each reaction family, the template was created based on known geometries of transition states from the literature or from manual calculations. The resulting Arrhenius parameters and rate coefficients at 1000K are given in Table 1.

These rate coefficients have been compared to experimental and theoretical values in literature and an excellent agreement is found, which is the same agreement found when calculating rate coefficients from CBS-QB3 calculations “manually”. Most of the rate coefficients are within a factor of 2 of experimental and theoretical literature rate coefficients.

Table 1: Arrhenius parameters and rate coefficients at 1000K of the eight reactions automatically calculated using Genesys. The first line of each reaction corresponds to the forward reactions, the second line to the reverse reaction. The pre-exponential factor and rate coefficients are expressed in s^-1 for monomolecular reactions and m³ mol^-1 s^-1 for bimolecular reactions. The activation energies are expressed in kJ mol^-1.

Conclusions

New methods have been developed to automatically calculate rate coefficients from ab initio calculations. For this, Genesys is used, which provides a connectivity based representation of the species, transition states and reactions. This information is used to create 3D coordinates of the atoms in the species and transition states. These coordinates are inputted to quantum chemistry calculations, yielding high level thermodynamic and high pressure limit kinetic data for gas-phase reactions.

The new methods are of great value to calculate a large number of rate coefficients from ab initio calculations without the need for manual interventions or user knowledge. This is especially of great interest for (1) building approximation methods for rate coefficient calculations such as group additivity models and (2) increasing the accuracy of important rate coefficients in a kinetic model.

                                                           

1.            Vandewiele NM, Van Geem KM, Reyniers MF, Marin GB. Genesys: Kinetic model construction using chemo-informatics. Chem Eng J. 2012;207:526-538.

2.            Steinbeck C, Han Y, Kuhn S, Horlacher O, Luttmann E, Willighagen E. The Chemistry Development Kit (CDK):  An Open-Source Java Library for Chemo- and Bioinformatics. Journal of Chemical Information and Computer Sciences. 2003;43(2):493-500.

3.            Halgren TA. Merck molecular force field .1. Basis, form, scope, parameterization, and performance of MMFF94. Journal of Computational Chemistry. 1996;17(5-6):490-519.

4.            Gaussian 09 [computer program]. Wallingford, CT, USA: Gaussian, Inc.; 2009.