# (571d) First Principles Based Automated Kinetic Model Generation Using on-the-Fly Ab Initio Calculations

#### AIChE Annual Meeting

#### 2017

#### 2017 Annual Meeting

#### Catalysis and Reaction Engineering Division

#### Reaction Engineering for Combustion and Pyrolysis

#### Wednesday, November 1, 2017 - 1:36pm to 1:58pm

**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^{1}, 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^{2}. 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)^{3}, 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^{4} 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^{3} 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.