Fostering automation in chemical kinetics: a protocol for bond energy computation and the implementation of a hierarchical approach for thermochemistry calculations

Automatic PES exploration and thermochemical parameters estimation are fundamental challenges to overcome to fully automate calculation procedures in predictive chemical kinetics from first principles. Our in-house software EStokTP was born with this final purpose and is continuously updated to foster automation and minimize human time spent on the calculations and thus the possibility of human mistakes. In the present work, three main protocols written in Python that work in synergy with EStokTP are presented: 1) the automatic generation of input data for EStokTP calculations, 2) the estimation of bond energy and 3) thermochemical parameters estimation in the form of NASA polynomials. The automatic generation of input data avoids the generation of the Z-matrix by the user, which is sometimes source of error caused by non-trivial geometry definition (bad first guess structure) or bad atom order in the Z-matrix (e.g., specific dihedral angle definition), and it is a tedious process if the required calculations involve a high number of species. The bond energies are calculated by generating the possible fragments obtained by the rupture of a bond of the original molecule; the information can be used for the selection of the main reaction channels to explore in a more detailed way and avoiding side reaction channels with too high activation energy. The algorithm was tested on 1,3-butadiene-2-ol and simulation results are compared with correspondent experimental evaluations. In particular, the estimated bond energy for fragmentation number 4 (H-abstraction from oxygen atom) is only 0.9 [kcal mol-1] higher than the latest literature theoretical estimate. Estimation of C_P^0 (T) and S^0 (T) is made by exploiting molecular translational, vibrational, and rotational partition function contributions and 1D hindered rotor theory. H^0 (T) is evaluated by estimating ∆H^0 (0 K) using the Connectivity Based Hierarchy method, which is a rung-based atomization scheme. ∆H^0 (0 K) is then corrected to ∆H^0 (298.15 K), with an extrapolation scheme based on experimental atomization enthalpies. Finally, a series of regressions evaluate the set of 14 coefficients required for the computation of NASA polynomials. The estimation of ∆H^0 (0 K) is tested on a set of 142 species, reporting a mean absolute error of 0.39 [kcal mol-1], while the influence of the level of theory and the correction of ∆H^0 (0 K) to ∆H^0 (298.15 K) are tested on a set of 8 molecules with a mean absolute error of 0.54 [kcal mol-1]. Finally, the NASA polynomials of isoprene and 1,3-butadiene-2-ol are evaluated and confronted with an external database result.

L’esplorazione automatica delle superfici di potenziale e la stima dei parametri termochimici sono sfide fondamentali da superare per automatizzare completamente le procedure di calcolo nella cinetica chimica predittiva. Il software EStokTP è nato con questo scopo finale ed è continuamente aggiornato per favorire l’automazione e ridurre al minimo il tempo speso dall’utente peri calcoli e quindi la possibilità di errori umani. Nel presente lavoro vengono presentati tre protocolli principali scritti in Python che lavorano in sinergia con EStokTP: 1) la generazione automatica di dati di input per i calcoli di EStokTP, 2) la stima dell'energia di legame e 3) la stima dei parametri termochimici sotto forma di polinomi NASA. La genereazione automatica degli input evita la scrittura a mano della matrice z da parte dell’utente, che talvolta risulta essere fonte di errori causati dalla definizione non banale della geometria (struttura di primo tentative errata) o dall’ordine di definizione degli atomi nella Z-matrix (per esempio nella definizione degli angoli diedri) e risulta essere un processo tedioso se i calcoli richiesti includono un elevato numero di specie dalla geometria complessa. Le energie di legame sono calcolate generando tutti i possibili frammenti ottenuti dalla rottura di un unico legame della molecola di partenza; le informazioni possono essere usate per la selezione dei principali canali di reazione da esplorare in maniera più dettagliata, evitando i canali di reazione secondari aventi una energia di attivazione troppo elevata. L’algoritmo è stato testato sul 1,3-butadiene-2-olo, e le simulazioni comparate con i risultati sperimentali in letteratura. In particolare, l’energia stimata per l’estrazione dell’idrogeno dal gruppo ossidrile è stata sovrastimata di sole 0.9 [kcal mol-1] rispetto alla più recente stima teorica in letteratura. La stima di C_P^0 (T) e S^0 (T) è effettuata utilizzando i contribute traslazionali, vibrazionali e rotazionali delle funzioni di partizione molecolare, e il trattamento di alcune vibrazioni come rotori impediti nella teoria dei rotori monodimensionali. H^0 (T) è valutata stimando dapprima ∆H^0 (0 K) usando un metodo basato sulla gerarchia di connettività, che è uno schema di atomizzazione a più livelli. Successivamente ∆H^0 (0 K) è corretto per valutare ∆H^0 (298.15 K), utilizzando uno schema di estrapolazione basato sulle entalpie sperimentali dei singoli atomi. Una serie di regressioni calcola infine il set di 14 coefficienti richiesti per i polinomi NASA. La stima del ∆H^0 (0 K) è stata testata su 142 specie, riportando un errore assoluto medio di 0.39 [kcal mol-1], mentre l’influenza del livello di teoria usato e la correzione da ∆H^0 (0 K) a ∆H^0 (298.15 K) è stata testata su un set di 8 specie, con un errore medio assoluto di 0.54 [kcal mol-1]. Infine i polinomi NASA di isoprene e 1,3-butadiene-2-olo sono stati valutati e confrontati con i risultati forniti da un database esterno.