Computational development of models and tools for the kinetic study of astrochemical gas-phase reactions

This PhD thesis focuses on the application and development of computational tools and methodologies for the modeling of the kinetics of gas-phase reactions of astrophysical interest in the interstellar medium (ISM). The complexity related to the investigation of chemical reactivity in space is mostly due to the extreme physical conditions of temperature, pressure and exposure to high-energy radiation, which in turn also lead to the formation of exotic species, like radicals and ions. Nevertheless, there is still much to be understood about the formation of molecules, the major issue being the lack of sufficient laboratory (experimental and computational) studies. A more detailed and accurate study of all the chemical processes occurring in the ISM will allow us to obtain the data necessary to simulate the chemical evolution of an interstellar cloud over time using kinetic models including thousands of reactions that involve hundreds of species. The collection of the kinetic parameters required for the relevant reactions has led to the growth of different astrochemical databases, such as KIDA and UMIST. However, the data gathered in these catalogues are incomplete, and rely extensively on crude estimations and extrapolations. These rates are of paramount importance to get a better comprehension of the relative abundances of the chemical compounds extrapolated by the astronomers from the spectral data recorded through the radio telescopes and the in-orbit devices, like the satellites. Accurate state-of-the-art computational approaches play a fundamental role in analyzing feasible reaction mechanisms and in accurately predicting the associated kinetics. Such approaches usually rely on chemical intuition where a by-hand search of the most likely pathways is performed. Unfortunately, thisprocedure can lead to overlook significant mechanisms, especially when large molecular systems are investigated. Increasing the size of a molecule can also increase the number of its possible conformers which can show a different chemical reactivity with respect to the same chemical partner. This brings to get very complex chemical reaction networks in which hundreds of chemical species are involved and thousands of chemical reactions can occur.During the last decades, a lot of effort has been done to develop computational techniques able to perform extensive and thorough investigations of complex reaction mechanisms. Such approaches rely on automated computational protocols which drastically decrease the risk of making blunders during the search for significant reaction pathways.Furthermore, the accurate characterization of the potential energy surfaces (PESs) critical points, like reactants, intermediates, transition states and products involved in the reaction mechanism, is crucial in order to carry out a reliable kinetic investigation. The kinetic analysis of an erroneous potential energy surface, would lead to gross errors in the estimation of the rate constants of the chemical species involved in the reaction.In order to avoid such errors, the combination of high-level electronic structure calculations via composite scheme can be helpful to get a more precise estimation of the energy barriers involved in the reaction mechanism. It has been proven that "cheap"[1] composite schemes can achieve subchemical accuracy without any empirical parameters and with convenient computation times, making them perfect for the purpose of this thesis.In recent decades, many efforts have been made to develop theoretical and computational methodologies to perform accurate numerical simulations of the kinetics of such complex reaction mechanisms in a wide range of thermodynamic conditions that mimic extreme reaction environmentsas for combustion systems, the atmosphere and the ISM. Such methodologies are based on the ab initio-transition-state-theory-based master equation approach, which allows the determination of rate coefficients and branching ratios of chemical species involved in complex chemical reactions. This methodology allows to make accurate predictions of the relative abundances of the reaction products for complex reactions even under conditions of temperature and pressure not experimentally accessible, such as those that characterize the ISM. Based on these premises, this dissertation has been focused on the application of a computational protocol for the ab initio-based computational modeling and kinetic investigation of gas-phase reactions which can occur in the ISM.This protocol is based on the application of validated methodologies for the automated discovery of complex reaction mechanisms by means of the AutoMeKin[2] program, the accurate calculation of the energetic of the potential energy surfaces (PESs) through the junChS and junChS-F12a "cheap" composite schemes and the kinetic investigation using the StarRate computer program specifically designed to study gas-phase reactions of astrochemical interest in conjunction with the MESS program. Furthermore, this dissertation has been also focused on the development and implementation of StarRate, a computer program for the accurate calculation of kinetics through a chemical master equation approach of multi-step chemical reactions. StarRate is an object-based program written in the so-called F language. It is structured in three main modules, namely molecules, steps and reactions, which extract the properties needed to calculate the kinetics for the single-step reactions partecipating in the overall reaction. Another module, in_out, handles program’s input and output operations. The main program,starrate, controls the sequences of the calling of the procedures contained in each of the three main modules.Through these modular structure, StarRate[3] can compute canonical and microcanonical rate coefficients taking into account for the tunneling effect and the energy-dependent and time-dependent evolution of the species concentrations involved in the reaction mechanism. Such protocol has been applied to investigate the formation reaction mechanisms of some complex interstellar polyatomic molecules, named interstellar complex organic molecules (iCOMs). More specifically, the formation of prebiotic iCOMs in space has raised considerable interest in the scientific community, because they are considered as precursors of more complex biological systems involved in the origin of life in the Universe. Debate on the origins of these biomolecular building blocks has been further stimulated by the discovery of nucleobases and amino acids in meteorites and other extraterrestrial sources. However, few insights on the chemistry which brings to the formation of such compounds is known.  References: [1] Jacopo Lupi,Silvia Alessandrini,Cristina Puzzarini,and Vincenzo Barone.junchs and junchs-F12 models:Parameter-free efficient yet accurate compositeschemes for energies and structures of noncovalent complexes. Journal of Chem-ical Theory and Computation, 17(11):6974–6992, 2021. PMID: 34677974.[2] Emilio Martínez-Núñez, George L. Barnes, David R. Glowacki, Sabine Kopec,Daniel Peláez, Aurelio Rodríguez, Roberto Rodríguez-Fernández, Robin J. Shan-non, James J. P. Stewart, Pablo G. Tahoces, and Saulo A. Vazquez.Au-tomekin2021: An open-source program for automated reaction discovery. Journalof Computational Chemistry, 42(28):2036–2048, 2021.[3] Surajit Nandi, Bernardo Ballotta, Sergio Rampino, and Vincenzo Barone.Ageneral user-friendly tool for kinetic calculations of multi-step reactions withinthe virtual multifrequency spectrometer project. Applied Sciences, 10(5), 2020

Gas-phase formation and isomerization reactions of cyanoacetaldehyde, a prebiotic molecule of astrochemical interest

Cyanoacetaldehyde (NC-CH2CH=O) is considered, together with guanidine and urea, as a precursor of the pyrimidine bases cytosine and uracil. Although it has not yet been detected in the interstellar medium (ISM), several hypotheses have been put forward about its synthesis in solution and in the gas phase. In this paper, we present a gas-phase model of the barrierless reaction between formyl (HCO) and cyanomethyl (CH2CN) radicals leading to cyanoacetaldehyde and focus on its evolution through isomerization and dissociation pathways. The potential-energy surface for all reactions has been explored by DFT calculations employing double-hybrid functionals and further refined through the "Cheap" composite scheme. Our results indicate that the direct association of the two reacting radicals (HCO and CH2CN) is strongly exothermic and thus thermodynamically favored under the harsh conditions of the ISM. Microcanonical rate constants computed with the help of the StarRate program for energies up to 6 kJ mol(-)(1) above the dissociation limit show that the most abundant products are the two conformers of cyanoacetaldehyde (nitrile and carbonyl groups in a cis or trans configuration) which, despite having comparable stability, are obtained with a cis/trans ratio of 0.35:0.65. The formation of other products with relative abundances not exceeding 10% is also discussed

Understanding fluorescent amyloid biomarkers by computational chemistry

Protein misfolding diseases, including neurodegenerative disorders like Alzheimer’s disease, are characterized by the involvement of amyloid aggregation, which emphasizes the need for molecular biomarkers for effective disease diagnosis. The thesis addresses two aspects of biomarker development: firstly, the computation of vibrationally resolved spectra of small fluorescent dyes to detect amyloid aggregation, and secondly, the binding and unbinding processes of a novel ligand to the target protein. In relation to the first aspect, a hybrid model for vibrational line shapes of optical spectra, called VCI-in-IMDHO, is introduced. This model enables the treatment of selected modes using highly accurate and anharmonic vibrational wave function methods while treating the remaining modes using the approximate IMDHO model. This model reduces the computational cost and allows for the calculation of emission line shapes of organic dyes with anharmonicity in both involved electronic states. The interaction between the dyes and their environment is also explored to predict the photophysical properties of the oxazine molecules in the condensed phase. The position and the choice of the solvent molecule have a significant impact on the spectra of the studied systems as they altered the spectral band shape. However, further studies are necessary to confirm the findings. In addition to neurodegenerative diseases, the systemic amyloidoses represent another group of disorders caused by misfolded or misassembled proteins. In the cardiac domain, the accumulation of amyloid fibrils formed by the transthyretin (TTR) protein leads to cardiac dysfunction and restrictive cardiomyopathy. The investigation of binding and unbinding pathways between the TTR protein and its ligands is crucial for gaining a comprehensive understanding and enabling early detection of systemic amyloidoses and related disorders. Hence, exploring the different binding modes and the dissociation pathways of TTR-ligand complex is the primary objective of the second aspect of this thesis. The experimental study provides evidence of binding and X-ray crystallographic structure data on TTR complex formation with the fluorescent salicylic acid-based pyrene amyloid ligand (Py1SA). However, the electron density from X-ray diffraction did not allow confident placement of Py1SA, possibly due to partial ligand occupancy. Molecular dynamics and umbrella sampling approaches were used to determine the preferred orientation of the Py1SA ligand in the binding pocket, with a distinct preference for the binding modes with the salicylic acid group pointing into the pocket.Deutsche Forschungs-gemeinschaft (DFG)/Emmy Noether/KO 5423/1- 1/E

Decomposition of Intermolecular Interactions in Ab Initio Spectroscopy

Spectroscopy, the molecular response to electromagnetic radiation of different wavelengths, is one of the most powerful experimental tools for interrogating a molecule's structure and dynamics as it interacts with its environment. However, relating a spectroscopic signature to a molecular picture relies on sophisticated computational approaches, which offer a wealth of methods for identifying structures, intermolecular interactions, and their correlation with spectroscopic response. This thesis focuses on the how to correlate a molecule's structure and interactions with its environment via ab initio calculation of spectroscopic parameters. To build a molecular picture of carbon dioxide dynamics in ionic liquids (ILs), quantum chemical calculations on small clusters qualitatively reproduced the experimental ordering for carbon dioxide's asymmetric vibrational stretch peak position which shifts when dissolved in a series of ILs with varying anions. To uncover the physical origin of the shift, the language of decomposition analysis based on absolutely localized molecular orbitals (ALMO-EDA) was translated from energies to vibrational frequencies. Geometric distortion of carbon dioxide, as a result of charge transfer (CT) from the anion into the carbon dioxide, is the driving force for differentiating the carbon dioxide asymmetric stretch shift in different IL anions. After validating these simple models, we further decomposed the CT contribution into geometry and curvature mechanisms, finding that CT is a significant contributor in both the geometry optimization and frequency calculation steps. A comparison between ALMO-EDA and symmetry-adapted perturbation theory (SAPT) showed that while dispersion dominates the binding energy, excellent correlation between both total interaction energies and individual components for ALMO-EDA and SAPT validates the use of DFT, enabling the construction of a semiempirical spectroscopic map. This decomposition presented the first application of an EDA outside the energy realm into molecular properties; however, it is not generally applicable to arbitrary perturbations. A reformulation of the canonical linear response equations for use with ALMOs provides a direct connection between EDA terms and their corresponding contribution to spectra. Results for argon-lithium cation dimer polarizabilities show that allowing CT is equally important in both the underlying ground-state wavefunction and the response calculation, and should not be confused with basis set superposition error

Simulations of proton transfer processes using reactive force fields

Within this thesis we presented the development of reactive force fields that are ca- pable to describe the dynamics of proton and hydrogen atom transfer processes. The presented implementation in CHARMM overcomes the limitation that bond break- ing and formation cannot be investigated by conventional classical MD simulations. Derived from high-level ab initio calculations this approach combines the accuracy of such calculations with the speed of MD simulations. The high-quality force fields of the prototype systems are comparable to high-level ab initio calculations in terms of structure and energy barriers. The PESs of proton transfer reactions are extremely sensitive with respect to the chemical environment. Nevertheless one is always able to classify the PT under investigation into symmetric and asymmetric PES. We de- veloped a series of parameter sets that are not only able to describe symmetric and asymmetric correctly but also can accommodate to different locations of energetic minima and barriers. The chosen three-dimensional potential energy functions have shown to be quite flexible and transferable in characterizing PT reactions in quite diverse chemical systems. The morphing transformation of MMPT force field param- eter, starting from one of our prototype systems to develop a new force field for a new molecular system which exhibits a similar topology in the PES along the proton reac- tion coordinates, has been shown to be successfully applicable in various examples. Energy scaling has been employed to investigate the effect on the proton transfer os- cillation in NH+ 4 · · ·NH3. New parameters through morphing have been developed for protonated diglyme, as well as for double proton transfer in 2PY2HP and for as- partic acid and water as model system for PT reactions in the active site of bacterial ferredoxin I. We applied the MMPT force field to investigate the vibrational infrared spectra of proton-bound species and explored the relationship of infrared spectra for protonated water dimer and protonated diglyme. The results for protonated water dimer compared well with other high-level calculations. Besides the further system- atic development of the morphing approach one can also employ the force field in combination with Feynman path integral methods. The MMPT force field could be a viable alternative to lower level quantum mechanical methods because the accuracy of the force field is only limited by the initial determination of the underlying PES for the PT of interest

Laboratory astrochemistry of dust and ice

Thin film growth and desorption behaviour of simple molecules on interstellar dust grain analogue surfaces has been investigated using a range of surface science techniques including temperature programmed desorption (TPD), reflection-absorption infrared (RAIR) and reflection-adsorption UV-Visible spectroscopy. The systems investigated use amorphous silica (aSiO2) as a mimic for bare interstellar dust grains and thin adsorbed films of ammonia (NH3), benzene (C6H6), carbon monoxide (CO), compact and porous amorphous solid water (c-ASW and p-ASW) crystalline solid water (CSW), methanol (CH3OH) and methyl formate (HCOOOCH3, MF). The optical properties for benzene (C6H6) were investigated using a newly designed and constructed UV/Visible spectrometer. Preliminary measurements of C6H6 on a highly-orientated pyrolytic graphite (HOPG) surface give the refractive index (n) as 1.43 ± 0.07 for a film of thickness (d) 261±5 nm. MF on aSiO2 was investigated using TPD, RAIRS and ab initio calculations. The TPD of MF is consistent with wetting of the aSiO2 surface. The binding energy of the monolayer was found to be 29.8±0.1 kJ mol-1 and that of the multilayer is 26.4±5.5 kJ mol-1 . This indicates that MF coupling to the aSiO2 surface is weak and only slightly stronger than the MF interaction with itself. Below 95 K, MF is in an amorphous phase and above 95 K, it is crystalline. A combination of measurements of spontaneous dipole orientation and RAIR spectra with computational chemistry supports the idea that the basis motif of the lattice in crystalline cis-MF is a ring dimer structure. A simple method was developed to synthesise the vibrational line profile of CO on a heterogeneous surface. The procedure developed allows the conversion of a distribution of binding energies, Edes, into a continuous distribution of vibrational frequencies, which can in turn be compared with experimental RAIRS data. The interaction of CO with a range of astrophysically surfaces including CH3OH, CSW, c-ASW, amorphous silica and NH3 on the aSiO2 substrate was investigated using TPD. Extended Inversion Analysis was used to determine the pre-exponential factor, distribution of Edes and the entropy of activation (∆ǂS) for desorption of CO from each surface

Theoretical Modeling of Water and Aqueous Systems

Thermodynamic and kinetic properties of liquid water and a variety of aqueous systems were investigated by means of theoretical methods. Focus was set on the prediction of mixing-induced changes in thermodynamic potentials, as well as on interfacial water dynamics. One of the employed methods is Quantum Cluster Equilibrium (QCE) theory, which introduces quantum chemistry to the well-established class of mixture theories. In this thesis, the origin of the two empirical QCE parameters was analyzed in great detail and their tight connection to the van der Waals equation of state was established. This connection was exploited to refine the binary mixture version of QCE theory. The proposed measures can eliminate the need for binary reference data and were proven to be viable approximations in a study on three different amide/water mixtures. Predicting miscibilities without need for binary reference data or resorting to empiricism is an important and open task in modern theoretical research on binary mixtures that is now within the reach of QCE theory. Like in any mixture theory, the success of QCE calculations depends sensitively on the choice of representative clusters, thus a systematic and empiricism-free scheme to generate cluster sets was proposed. Therein, only global minimum structures, obtained from a genetic-algorithm based global geometry optimization, are employed. Such cluster sets reach almost chemical accuracy in the prediction of excess enthalpies of mixing. However, correctly describing excess entropies turned out to be out of their reach, due to the lack of energetically less stable, yet entropically important cluster structures. Suitable extensions of the scheme addressing this issue have been outlined. QCE theory can go beyond the prediction of thermodynamic potentials and other structural insight into liquids, as well. In this thesis, the ionic product of water was calculated by QCE theory, and thus the path to the investigation of various other acid-base related phenomena has been opened. The proposed theoretical improvements and the performed QCE calculations were made possible by the Peacemaker QCE software, which was rewritten from scratch as part of this work. The design, implementation, and use of the code are documented herein. Some of the most fundamental properties of water and aqueous systems are dynamic in their nature, and thus beyond the reach of QCE theory. In the spirit of modern multimethod research, atomistic simulation was the second method of choice employed in this thesis and was applied to a selection of problems that occur on the nanoscale. A large part of this work was devoted to hydrogen bond and allied dynamics, which are a major driving force for processes occurring in liquid water. Further focus was set on electric field effects, which influence various applications ranging from nanofluidic devices to membrane ion channels. On the nanoscale, the spontaneous orientational polarization of water can couple with electric field alignment, resulting in an asymmetric behavior at opposing surfaces—a situation that has previously been described as field-induced Janus interface. Here, a new and significant field polarity (sign) dependence of the dipolar reorientation dynamics in water hydration layers was uncovered. Imposition of an electric field across a nanopore can lead to differences in response times of interfacial water polarization of up to two orders of magnitude, with typical time scales being in the picosecond regime. Coupling between interfacial polarization and interfacial density relaxations was revealed, as well. The surprisingly strong asymmetry in the dynamic response at opposing surfaces is even more pronounced than the known static properties of a field-induced Janus interface. Cavities found in nature and technology are often spherical, and water dynamics in such nanoconfinement was investigated, as well. A diffusive model was constructed by Bayesian inference from simulation data, which describes the single-particle dynamics of water molecules inside spherical cavities (fullerenes). The propagators of the diffusion model show good agreement with simulation data over four orders of magnitude, instilling great confidence in the model. There was no a priori reason to believe in the existence of such a diffusion model, but after having established its validity, hydrogen bond kinetics could be meaningfully treated within the same diffusion model that applies to bulk water. Overall, hydrogen bond lifetimes slow down with decreasing cavity size. An attempt was made to predict hydrogen bond time correlation functions from a simple pair diffusion equation with sink and source terms corresponding to hydrogen bond breaking and formation, but the model could not be found to be reliable in spherical nanoconfinement. Various ways to improve upon this procedure have been proposed for follow-up studies