(571d) 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.