International Journal of Molecular Sciences Article Gas-Phase vs. Grain-Surface Formation of Interstellar Complex Organic Molecules: A Comprehensive Quantum-Chemical Study Berta Martínez-Bachs and Albert Rimola * Departament de Química, Universitat Autònoma de Barcelona, 08193 Bellaterra, Catalonia, Spain; berta.martinez@uab.cat * Correspondence: albert.rimola@uab.cat; Tel.: +34-93-581-3723 Abstract: Several organic chemical compounds (the so-called interstellar complex organic molecules, iCOMs) have been identified in the interstellar medium (ISM). Examples of iCOMs are formamide (HCONH2 ), acetaldehyde (CH3 CHO), methyl formate (CH3 OCHO), or formic acid (HCOOH). iCOMs can serve as precursors of other organic molecules of enhanced complexity, and hence they are key species in chemical evolution in the ISM. The formation of iCOMs is still a subject of a vivid debate, in which gas-phase or grain-surface syntheses have been postulated. In this study, we investigate the grain-surface-formation pathways for the four above-mentioned iCOMs by transferring their primary gas-phase synthetic routes onto water ice surfaces. Our objective is twofold: (i) to identify potential grain-surface-reaction mechanisms leading to the formation of these iCOMs, and (ii) to decipher either parallelisms or disparities between the gas-phase and the grain-surface reactions. Results obtained indicate that the presence of the icy surface modifies the energetic features of the reactions compared to the gas-phase scenario, by increasing some of the energy barriers. Therefore, the investigated gas-phase mechanisms seem unlikely to occur on the icy grains, highlighting the distinctiveness between the gas-phase and the grain-surface chemistry. Keywords: astrochemistry; interstellar grains; water ice surfaces; molecule formation; iCOMs; computational chemistry Citation: Martínez-Bachs, B.; Rimola, A. Gas-Phase vs. Grain-Surface 1. Introduction Formation of Interstellar Complex The interstellar medium (ISM) is the region between stars. Although not visible to the Organic Molecules: A Comprehensive naked eye, ISM contains matter in the form of gas-phase molecules and solid-state dust Quantum-Chemical Study. Int. J. Mol. Sci. 2023, 24, 16824. https://doi.org/ grain particles (with a mass proportion of 99% and 1%, respectively). In dense molecular 10.3390/ijms242316824 clouds (regions where star formation begins and is featured to have temperatures of 5–10 K and gas density of 104 cm−3 ), interstellar grains are sub-micrometer-sized particles, Academic Editor: Dongho Kim constituted by a refractory core (of silicates or carbonaceous materials) that is covered by Received: 27 October 2023 ice mantles, predominantly of H2 O but also of other volatile species like CO, CO2 , NH3 , Revised: 20 November 2023 and CH3 OH, all of them identified using near-infrared observational measurements [1,2]. Accepted: 24 November 2023 Regarding the gas-phase components, to date, approximately 300 interstellar species have Published: 27 November 2023 been identified via rotational emission spectroscopy in the millimeter and sub-millimeter range [3]. Among these species, ca. one-third belongs to the group called interstellar complex organic molecules (iCOMs), which are carbon-bearing species containing at least five atoms [4–6]. iCOMs can be precursors of more-complex organic species, in which some Copyright: © 2023 by the authors. of them can be of biological interest [7–10]. Therefore, the study of iCOMs formation is Licensee MDPI, Basel, Switzerland. important not only for understanding the chemical complexity of space but also for offering This article is an open access article valuable insights on the potential origins of organic compounds that could be relevant to distributed under the terms and the chemistry of life. conditions of the Creative Commons Currently, there are two prevailing paradigms for explaining the formation of iCOMs Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ in the ISM: either through gas-phase reactions or on the surface of interstellar grains [11–16]. 4.0/). In both theories, the first step advocates the presence of hydrogenated species in the icy Int. J. Mol. Sci. 2023, 24, 16824. https://doi.org/10.3390/ijms242316824 https://www.mdpi.com/journal/ijms Int. J. Mol. Sci. 2023, 24, 16824 2 of 25 grains (e.g., H2 O, NH3 , CH4 , H2 CO, or CH3 OH), which are formed by hydrogenation (i.e., H addition) onto atoms and simple species present on the grains [17]. For instance, H addition to O, N, and CO can give rise to the formation of H2 O, NH3, and CH3 OH, respectively [18–20]. In the gas-phase scenario, the frozen hydrogenated species partially desorb into a gas-phase through thermal or non-thermal processes, where they react with other gaseous species forming iCOMs [13,14,21,22]. Due to the extremely low temperatures, feasible pathways in these harsh environments require both exothermic reactions and barrierless mechanisms or, at least, pathways presenting low energy barriers. Therefore, gas-phase reactions involve either ions or radicals, and can be categorized into two main groups: ion-neutral or neutral–neutral reactions. Regarding the neutral–neutral reactions, they are reactions between a radical and a closed-shell species, and those that are not barrierless are characterized to form a pre-reactive complex [16]. At very low temperatures, the prereactive complex can remain stable long enough to allow the system to overcome the low energy barrier of the reaction [23]. Neutral–neutral reactions primarily lead to the formation of a radical intermediate (because of the reactive coupling of the two species), which, under strict gas-phase conditions, cannot dissipate its internal energy, and accordingly cannot be stabilized. Due to its high internal energy, this radical intermediate can either dissociate back or evolve, yielding newly formed species [14,21,22,24,25]. In the grain-surface scenario, the hydrogenated icy species trapped on the grain surfaces are photo-dissociated by ultra-violet (UV) photons generating radicals. These radicals become more mobile when the grain temperature increases to about 20–30 K during the collapse phase, allowing them to diffuse over the surface and recombine to form iCOMs, in which such a radical–radical coupling is assumed to be barrierless [11,12]. Finally, iCOMs are released to the gas phase using thermal or non-thermal processes. However, calculations have shown that radical–radical reactions on water ice surfaces can present barriers and can have competitive channels, such as H-abstraction, significantly limiting the efficiency of the processes [26–30]. Therefore, alternative processes towards iCOMs synthesis have started to be considered, such as non-diffusive processes [31,32] involving gas–grain reactions between radicals (gas) and components of the ices (grain) [27,33–35]. The mechanisms and pathways that explain the formation of iCOMs still constitute a vivid debate. However, it seems that both paradigms, the gas-phase and the grain-surface mechanisms, are essential to account for the diversity and abundances of iCOMs observed in the ISM [11]. Determining whether the formation of a particular iCOM is governed by the gas-phase or the grain-surface chemistry is a complex matter since it depends on the intrinsic nature of the target iCOM, and accordingly a case-by-case investigation is required [16]. In this work, we investigate, by means of quantum-chemical simulations, the forma- tion pathways of four astrochemically relevant iCOMs that are recurrently identified in different star-forming regions: formamide (HCONH2 ), acetaldehyde (CH3 CHO), methyl formate (HCOCH3 ), and formic acid (HCOOH). Our approach involves a comprehensive investigation of the formation routes of these species on interstellar water ice models. How- ever, the mechanisms under study in this contribution are inspired by the corresponding neutral–neutral gas-phase synthetic routes, which are considered to be major sources of these species in the gas phase. These gas-phase formation routes have not been considered on the ice mantles, and intriguingly, the interaction of the reactive species with the surfaces can modify the energetic features of the reactions, either increasing or decreasing the en- ergy barriers, and favoring or disfavoring their thermodynamics. Accordingly, our goal is twofold: (i) to assess whether these gas-phase formation routes can also be potential paths for their synthesis on interstellar water ice mantles, and (ii) to determine any analogy or lack thereof of the gas-phase mechanisms compared to reactions occurring on the surfaces of interstellar water ice mantles. Int. J. Mol. Sci. 2023, 24, 16824 3 of 25 2. Results 2.1. Gas-Phase Reactions Initially, the reactions were simulated under gas-phase conditions (namely, in the absence of the water ice models) using the CCSD(T)//DFT methodology (see Section 4.2). These primary calculations, in addition to allowing us to understand the intrinsic electronic properties of the reactions, also serve to carry out a benchmarking study to select the most accurate DFT functional that best describes each of them. Details of the benchmarking study can be found in Appendix A. It is worth noting that in the energetics of these gas-phase reactions, the ZPE-corrections added to the CCSD(T) potential-energy values are those computed at the DFT/6-311+G(d,p) theory level. 2.1.1. Formamide Formation Formamide (HCONH2 ) is an iCOM that was first detected in the ISM in 1971 [36]. Since then, it has been identified in several star-forming regions, becoming a species commonly found in such regions [37–43]. Formamide is an iCOM of great interest because it is the simplest iCOM containing the four most important elements from a biological point of view: H, C, N, and O. Furthermore, formamide is the simplest molecule containing an amide bond (-CO-NH-), which is the same bond that joins amino acids in peptides. Therefore, formamide represents a turn out in chemical complexity within the ISM, increasing its interest from the standpoint of prebiotic chemistry. Several pathways relative to the formation of formamide on grain surfaces have been proposed [27,29,44,45], but recently, a neutral–neutral gas-phase formation route for formamide has been postulated by Barone et al. to be particularly favorable [24]. It consists of a two-step reaction: the first step involving the coupling of formaldehyde (H2 CO) and the NH2 radical, leading to the formation of a radical intermediate (H2 CONH2 ), and the second step involving the dissociation of the radical intermediate yielding formamide and a hydrogen atom. That is, H2 CO + NH2 → H2 CONH2 → HCONH2 + H (1) This mechanism, although the PES presents potential energy barriers, was considered to be barrierless because most of the transition states are lower in energy than the initial reactants. In this mechanism, the radical H2 CONH2 intermediate was found to be the most stable species on the PES, and, thus, the products were less stable than the intermediate. However, after adding the ZPE corrections, the intermediate and the product exhibit similar stability. In some studies, the product is the most stable structure [24], while in others the intermediate is the most stable one [25,46]. Since this reaction is associated with the release of an H atom, kinetic calculations based on the Rice–Ramsperger–Kassel–Marcus (RRKM) theory demonstrated its feasibility in the ISM [25,46]. The formation of formamide following this Reaction (1) is characterized here as CCSD(T)//DFT, in which the M06-2X functional is found as the most accurate DFT method for this reaction (see Appendix A for details). Figure 1 displays the geometries and the ZPE-corrected PES for this reaction. The results obtained agree with the previous studies, identifying the intermediate radical I as the most stable structure of the pure PES. However, similar to previous studies [24,25,46], after adding the ZPE-corrections, the intermediate and the products exhibit similar stability, the product being 3.7 kJ mol−1 more stable than the intermediate. The first step of the reaction (the coupling between the reactants to form I) exhibits an intrinsic energy barrier (energy difference between TS1 and PR, see Figure 2) of 22.6 kJ mol−1 , in which the TS1 is 18.8 kJ mol−1 higher than the reactants R asymptote. The second step (intermediate dissociation) presents an intrinsic energy barrier (here, as energy difference between TS2 and I) of 45.6 kJ mol−1 , in which the TS2 is 11.2 kJ mol−1 higher than the reactants R asymptote. barrier (here, as energy difference between TS2 and I) of 45.6 kJ mol−1, in which the TS2 is Int. J. Mol. Sci. 2023, 24, 16824 11.2 kJ mol−1 higher than the reactants R asymptote. 4 of 25 Figure 1. CCSD(T)//M06-2X-D3 ZPE-corrected PES and relative energies (with respect to the R Figure 1. CCSD(T)//M06-2X-D3 ZPE-corrected PES and relative energies (with respect to the R as- asymptote), and M06-2X-D3-optimized geometries for the gas-phase formamide formation. Distances ymptote), and M06-2X-D3-optimized geometries for the gas-phase formamide formation. Distances are in Å and relative potential energies in kJ mol−1 . Color coding: white, H atoms; grey, C atoms; are in Å and relative potential energies in kJ mol−1. Color coding: white, H atoms; grey, C atoms; blue,blue, N atoms; and red, O atoms. N atoms; and red, O atoms. 2.1.2. Acetaldehyde Formation 2.1.2. Acetaldehyde Formation Acetaldehyde (CH3 CHO) is one of the most common iCOMs identified in the ISM, Acetaldehyde and first detected(CH 3CHO) is one of the most common iCOMs identified in the ISM, in 1971 [47]. Furthermore, acetaldehyde is a possible precursor of several andcarbohydrates, first detected thus in 1971 being[47]. Furthermore, an iCOM acetaldehyde of great prebiotic is a[48,49]. potential possible precursor Different of sev- acetalde- eralhyde carbohydrates, thus being formation routes an iCOM have been of great discussed prebiotic in the potential literature, [48,49]. including both Different gas-phase ac- etaldehyde synthesisformation and reactions routes have on occurring been the discussed surfaces of in the ice water literature, including both gas- mantles [26,28,35,45,50–54]. Skouteris phase et al. [21] synthesis andproposed a neutral–neutral reactions occurring on gas-phase mechanism the surfaces ofbased wateron the iceethanol mantles tree, in which ethanol [26,28,35,45,50–54]. (CH3 CH Skouteris 2 OH) et al. is proposed [21] consideredathe parent molecule neutral–neutral of different gas-phase iCOMs mechanism such based onastheglycolaldehyde ethanol tree, (HCOCH in which2 OH) formic ethanol (CHacid (HCOOH), acetic acid (CH2 COOH), 3CH2OH) is considered the parent mole- and also acetaldehyde. cule of different iCOMs such as glycolaldehyde (HCOCH2OH) formic acid (HCOOH), The formation of acetaldehyde from the ethanol tree mechanism involves two sequen- acetic acid (CH2COOH), and also acetaldehyde. tial reactions. The first reaction is the α H-abstraction from ethanol using an OH radical, The formation of acetaldehyde from the ethanol tree mechanism involves two se- resulting in the formation of the CH3 CHOH radical and water, quential reactions. The first reaction is the α H-abstraction from ethanol using an OH rad- ical, resulting in the formation of 2the CH3 CH OHCH 3CHOH + OH → CHradical 3 CHOHand + Hwater, 2O (2) CH3CH The H-abstraction reaction 2OH can + OH occur → CH at three 3CHOH different + Hresulting sites, 2O in the formation (2) of the radicals CH2 CH3 OH, CH3 CHOH, or CH3 CH3 O, and water. Theoretical studies inves- The H-abstraction tigating reaction this reaction, either can occur identified at three different the α H-abstraction transitionsites, stateresulting in the slightly below thefor- mation energyof the levelradicals CH2CH(2.5 of the reactants 3OH, kJ CH mol3−CHOH, 1 ) [55] oror CH3CH slightly 3O, and above water. (0.7 to 1.8 kJTheoretical stud- mol−1 ) [56,57]. ies investigating this reaction, either identified the α H-abstraction Additionally, several experimental investigations, conducted at temperatures relevant totransition state slightly below the energy atmospheric or level of the conditions, combustion reactants (2.5 have kJ concluded mol−1) [55]that or slightly above (0.7 to the α H-abstraction 1.8 kJ is the molprimary −1 channel ) [56,57]. for this reaction Additionally, several [58,59]. experimental investigations, conducted at tempera- The second tures relevant reaction is the to atmospheric or addition combustion of anconditions, oxygen atom haveto the CH3 CHOH concluded thatradical, the α H- resulting in the formation of a radical intermediate abstraction is the primary channel for this reaction [58,59]. CH 3 CHOOH, which subsequently dissociates The second to yield acetaldehyde reaction and anof is the addition OH anradical, oxygen atom to the CH3CHOH radical, re- sulting in the formation CH3of a radical CHOH +O→ intermediate CH3 CHOOHCH CHOOH, →3CH which subsequently(3)dis- 3 CHO + OH sociates to yield acetaldehyde and an OH radical, This second reaction was theoretically investigated by Skouteris et al. [21]. The first CH3CHOH + O → CH3CHOOH → CH3CHO + OH step of this second reaction was found to be truly barrierless, in which the coupling of the (3) 3 atomic O( P) to the CH3 CHOH radical occurred in a spontaneous way, forming directly This second reaction was theoretically investigated by Skouteris et al. [21]. The first the CH CHOOH intermediate. The second step (intermediate dissociation) presented step of this3 second reaction was found to be truly barrierless, in which the coupling of the an energy barrier but submerged below the energy of the reactants. The CH3 CHOOH atomic O(3P) to the intermediate was CH 3CHOH radical occurred in a spontaneous way, forming directly found to be the most stable structure of the PES, giving rise to overall the unfavorable CH3CHOOH intermediate. reaction TheAdditionally, energetics. second stepthe(intermediate dissociation) branching ratio presented an for the acetaldehyde energy barrier but submerged below the energy of the reactants. The CH CHOOH formation compared to other possible channels (i.e., formic acid and acetic acid formation), 3 inter- mediate was found was found to be0.075. to be only the most stable However, structure Vazart of the et al. [22] byPES, meansgiving rise to overall of RRKM-based unfa- kinetic calculations vorable reactionestimated the overall energetics. rate constant Additionally, forbranching the the acetaldehyde ratio formation, suggesting for the acetaldehyde Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 5 of 26 formation compared to other possible channels (i.e., formic acid and acetic acid for- Int. J. Mol. Sci. 2023, 24, 16824 5 of 25 mation), was found to be only 0.075. However, Vazart et al. [22] by means of RRKM-based kinetic calculations estimated the overall rate constant for the acetaldehyde formation, suggesting that this ethanol-tree-based path is kinetically favorable as a potential gas- that this phase ethanol-tree-based source for the formation pathofisinterstellar kinetically favorable as a potential gas-phase source for acetaldehyde. the formation of interstellar acetaldehyde. Here, we calculate this gas-phase formation route for acetaldehyde at the Here, we calculate this gas-phase formation route for acetaldehyde at the CCSD(T)//M06- CCSD(T)//M06-250 2X-D3 level (the most accurate methodology among the tested ones, 250 2X-D3 level (the most accurate methodology among the tested ones, see Appendix A see Appendix A for details). The geometries and the ZPE-corrected PESs are illustrated in for details). The geometries and the ZPE-corrected PESs are illustrated in Figure 2. For Figure 2. For reaction 2, the results identify the products P1 as the most stable structure of Reaction (2), the results identify the products P1 as the most stable structure of the profile. the profile. Although this Although reactionthis reaction exhibits exhibitsenergy an intrinsic an intrinsic barrierenergy barrier of 18.1 kJ mol−1of 18.1 kJdifference (energy mol−1 (en- ergy difference between TS1 and between PR), theTS1 and PR), transition theTS1 state transition state is slightly TS1in below is energy slightlythan belowthe in energy reactant than the reactant R1 asymptote − 1 (−1.3 kJ mol −1). In contrast, for reaction 3, the intermediate R1 asymptote (−1.3 kJ mol ). In contrast, for Reaction (3), the intermediate I is the most I stable is the most stable structure structure along along the the pathway, to pathway, the detrimentto theofdetriment of the final the final product P2. product The O(3 P)P2. The O( P)step addition3 addition step is barrierless, is barrierless, as it is a radical–radical as it is a radical–radical coupling, whilecoupling, while thestep the dissociation dissoci- of I ation stepan presents ofintrinsic I presents an intrinsic energy barrier of energy mol−1 of 72.8 kJbarrier 72.8 kJ (energy mol−1 (energy difference between difference be- TS2 and I). tween However,TS2 itand I). However,that is noteworthy it isthe noteworthy energy level that the transition of the energy levelstateofTS2 theistransition state significantly TS2 loweris significantly than that of thelower than that reactants of the reactants R2 asymptote. R2 asymptote. CCSD(T)//M06-2X-D3 Figure2.2.CCSD(T)//M06-2X-D3 Figure ZPE-corrected ZPE-corrected PESPES andand relative relative energies energies (with(with respect respect tocorre- to the the corresponding sponding R1 andR1R2and R2 asymptotes), asymptotes), and M06-2X-D3-optimized and M06-2X-D3-optimized geometries geometries for theforgas-phase the gas-phase acetal- acetaldehyde formation. Distances −1 . Color coding: dehyde formation. Distances are in are in Årelative Å and and relative potential potential energies energies in kJ in mol kJ−1mol . Color coding: white, Hwhite, atoms;H grey, C grey, atoms; atoms; and red, C atoms; andO red, atoms. O atoms. 2.1.3. Methyl Formate Formation Methyl formate (HCOOCH3 ), the simplest ester identified in the ISM, was first detected in 1975 [60]. Moreover, methyl formate is detected in more abundance than its structural isomers acetic acid (CH3 COOH) and glycolaldehyde (HCOCH2 OH) [61]. Methyl formate (HCOOCH3), the simplest ester identified in the ISM, was first tected in 1975 [60]. Moreover, methyl formate is detected in more abundance than its st tural isomers acetic acid (CH3COOH) and glycolaldehyde (HCOCH2OH) [61]. Int. J. Mol. Sci. 2023, 24, 16824 Various mechanisms for the formation of methyl formate have been studied 6 of 25 in re years [11,62,63]. A neutral–neutral gas-phase synthetic pathway was proposed by Ru [64] and Balucani [14], pointing out dimethyl ether (CH3OCH3) as a gas-phase precu for methyl Various formate. mechanisms The forproposed mechanism, the formation previously of methyl formate haveinvestigated been studied in byre-Hoyerm [65] cent and yearsSong [66], involves [11,62,63]. a two-step A neutral–neutral reaction. gas-phase The first synthetic step iswas pathway theproposed additionby of an oxy Ruaud [64] and Balucani [14], pointing out dimethyl ether (CH OCH ) as a gas-phase atom to the CH3OCH2 radical, which is barrierless due to its radical–radical coupling 3 3 precursor ture, for methyl resulting in theformate. The proposed formation mechanism, of H2COOCH previously investigated by Hoy- 3 as a radical intermediate. The second ermann [65] and Song [66], involves a two-step reaction. The first step is the addition of consists of the dissociation of this radical intermediate yielding methyl formate and a an oxygen atom to the CH3 OCH2 radical, which is barrierless due to its radical–radical drogen couplingatom, nature, resulting in the formation of H2 COOCH3 as a radical intermediate. The second step consists of the dissociation of this radical intermediate yielding methyl formate CH2OCH3 + O → H2COOCH3 → HCOOCH3 + H and a hydrogen atom, Dedicated calculations CH2 OCH3 + found O → H2the radical intermediate as the most stable (4) specie COOCH 3 → HCOOCH3 + H the path. Additionally, different competitive channels were explored, including the di Dedicated ciation calculations of the radical found theinto intermediate radical intermediate several productsassuch the most as CHstable species 3O and CH2O, CH2 of the path. Additionally, different competitive channels were and CH2O, CH3 and HCOOH, HCO and CH3OH, and H2 and CH3CO. However, explored, including the kin dissociation of the radical intermediate into several products such as CH3 O and CH2 O, calculations indicated that the formation of methyl formate dominates over the other CH2 OH and CH2 O, CH3 and HCOOH, HCO and CH3 OH, and H2 and CH3 CO. However, tential kinetic competitive channels, calculations indicated inthe that which tunneling formation effects of methyl were underlined formate dominates oversince thethe reac involves the release other potential of a channels, competitive hydrogen in atom which [66]. tunneling effects were underlined since the Theinvolves reaction synthesis of methyl the release of a formate hydrogenfollowing atom [66]. this neutral–neutral gas-phase mechan Theinvestigated, is here synthesis of methyl in thisformate following case using this neutral–neutral gas-phase the CCSD(T)//B3LYP-D3(BJ) mecha- methodology (for nism is here investigated, in this case using the CCSD(T)//B3LYP-D3(BJ) methodology reaction, the B3LYP-D3(BJ) method was found to be the most accurate, see Append (for this reaction, the B3LYP-D3(BJ) method was found to be the most accurate, see Ap- for details). Figure 3 shows the optimized geometries and the corresponding ZPE- pendix A for details). Figure 3 shows the optimized geometries and the corresponding rected PES. PES. ZPE-corrected Figure 3. CCSD(T)//B3LYP-D3(BJ) ZPE-corrected PES and relative energies (with respect to the R Figure 3. CCSD(T)//B3LYP-D3(BJ) ZPE-corrected PES and relative energies (with respect to t asymptote), and B3LYP-D3(BJ)-optimized geometries for the gas-phase methyl formate formation. asymptote), and B3LYP-D3(BJ)-optimized geometries for the gas-phase methyl formate forma Distances are in Å and relative potential energies in kJ mol−1 . Color coding: white, H atoms; grey, Distances are in Å and C atoms; and red, O atoms. relative potential energies in kJ mol−1. Color coding: white, H atoms; gre atoms; and red, O atoms. Our results are consistent with the previous studies. The radical intermediate is the most Our stableresults structure of the are energy profile consistent with(more than the products the previous studies.P),The the atomic radicalO intermediate addition is step is barrierless, most stable structure and the intermediate of the energy dissociation profile step (more presents than thean intrinsic products energy P), barrier the atomic O a of 49.0 kJ mol−1 (energy difference between TS and I). Nevertheless, the transition state TS tion step is barrierless, and the intermediate dissociation step presents an intrinsic ene is significantly more submerged in energy than the reactants R asymptote. barrier of 49.0 kJ mol−1 (energy difference between TS and I). Nevertheless, the transi state 2.1.4. TS is significantly Formic more submerged in energy than the reactants R asymptote. Acid Formation Formic acid (HCOOH) is an organic molecule detected for the first time in the ISM 2.1.4. in 1971Formic Acidacid [67]. Formic Formation is the simplest carboxylic acid and serves as a fundamental com- pound for the formation of more-complex carboxylic acids, including biomolecules [68]. Formic acid (HCOOH) is an organic molecule detected for the first time in the ISM Several pathways for the formation of formic acid have been considered in the litera- 1971 [67]. Formic acid ture [11,14,64,69,70]. is the simplest The UMIST carboxylic [71] and KIDA acid and [72] databases serves report as a methods various fundamental c pound for the formation of more-complex carboxylic acids, including biomolecules of gas-phase synthesis for formic acid. One of them is the neutral–neutral gas-phase for- [ Several pathways for the formation of formic acid have been considered in the literatu [11,14,64,69,70]. The UMIST [71] and KIDA [72] databases report various methods of ga phase synthesis for formic acid. One of them is the neutral–neutral gas-phase formatio route, which is analogous to the two-step formamide formation route considered in th Int. J. Mol. Sci. 2023, 24, 16824 work. The first step involves the coupling between formaldehyde (H2CO) and a OH rad 7 of 25 cal, leading to the formation of a radical intermediate, H2COOH. The second step consis of the dissociation of the radical intermediate, forming formic acid and releasing a hydr mation gen route, which is analogous to the two-step formamide formation route considered atom, in this work. The first step involves the coupling between formaldehyde (H2 CO) and a H2CO + OH OH radical, leading to the formation → H2COOH of a radical → HCOOH intermediate, + HThe second step H2 COOH. ( consists of the dissociation of the radical intermediate, forming formic acid and releasing a Yetter atom, hydrogen [73] conducted experimentally studies on this reaction, although not und ISM conditions. They H found OH → 2 CO +that at H 2 COOH → HCOOH temperatures above+200 H K, a competitive(5)channel this mechanism, the H-abstraction Yetter [73] conducted leading experimentally toon studies HCO and H2O, this reaction, was the although notmain underbranching ISM r tio. D’Anna They conditions. [74] performed found that atatemperatures computational abovework 200 K,ona this reaction competitive obtaining channel of thisan energ barrier for the first transition state of 27.8 kJ mol2 ; although, they did notratio. mechanism, the H-abstraction leading to HCO and H O, −1 was the main branching consider t D’Anna [74] performed a computational formation of the pre-reactant complex. work on this reaction obtaining an energy barrier Finally, Xu [75] focused on the simulation of t for the first transition state of 27.8 kJ mol−1 ; although, they did not consider the formation first step of this reaction, finding an energy barrier of 23.8 kJ mol−1. Despite these mix of the pre-reactant complex. Finally, Xu [75] focused on the simulation of the first step results, of this different models reaction, finding anconsider this reaction energy barrier of 23.8 kJasmolone−1 . of the major Despite these channels towards fo mixed results, mic acid formation different [11,12,76]. models consider this reaction as one of the major channels towards formic acid Here, [11,12,76]. formation we explore this formation pathway at the CCSD(T)//B3LYP-D3(BJ) level, whichHere, we reaction for this explore this formation pathway B3LYP-D3(BJ) at the to was found CCSD(T)//B3LYP-D3(BJ) be the most accurate level, method in (see A which for this reaction B3LYP-D3(BJ) was found to pendix A for details). The obtained results are shown in Figure 4.be the most accurate method (see Appendix A for details). The obtained results are shown in Figure 4. Figure 4. CCSD(T)//B3LYP-D3(BJ) ZPE-corrected PES and relative energies (with respect to the Figure 4. CCSD(T)//B3LYP-D3(BJ) ZPE-corrected PES and relative energies (with respect to the R asymptote), and B3LYP-D3(BJ)-optimized geometries for the gas-phase formic acid formation. asymptote), and B3LYP-D3(BJ)-optimized geometries − for the gas-phase formic acid formation. D Distances are in Å and relative potential energies in kJ mol 1 . Color coding: white, H atoms; grey, tances are in Å and relative potential energies in kJ mol−1. Color coding: white, H atoms; grey, C atoms; and red, O atoms. atoms; and red, O atoms. The results are in very good agreement with the previous works. The most stable The results structure are inprofile of the energy very good agreement is the radical with the intermediate, previous I, which works. The is significantly most stab lower in energyofthan structure thethe final products energy profile isP.the Theradical couplingintermediate, step of the reaction presents I, which an intrinsic lower is significantly potential energy barrier (energy difference between TS1 and PR) of 33.1 kJ mol −1 (very energy than the final products P. The coupling step of the reaction−1presents an intrins similar to that reported in previous works), in which the TS1 is 18.7 kJ mol higher than potential energy barrier (energy difference between TS1 and PR) of 33.1 kJ mol−1 (ve the reactants R asymptote. The dissociation step presents an intrinsic energy barrier of similar 59.2 kJto that mol reported −1 (energy in previous difference betweenworks), TS2 andinI), which the but the TS2TS1 is 18.7lays structure kJ mol below higher −1 in th theenergy reactants R asymptote. with respect The dissociation to the reactants R asymptote. step presents an intrinsic energy barrier 59.2 kJ mol−1 (energy difference between TS2 and I), but the TS2 structure lays below 2.2. Reactivity on the Water Ice Surface Models energy with respect to the reactants R asymptote. Aiming to determine if the previously studied gas-phase reactions can also occur on the surfaces of water ice mantles, we computed them on our CWS and ASW periodic water 2.2. Reactivity on the Water Ice Surface Models ice models. It is worth mentioning that on the CWS, due to its inherent crystallinity, only oneAiming to determine surface reaction if the previously site is available. However,studied gas-phase this is not reactions the case for can the ASW also occur o model, theinsurfaces of water which different ice mantles, adsorption siteswe cancomputed them be identified, on our which CWS andreaction are potential ASW periodic sites. wat Nevertheless, computing all the reactions on all these surface sites is computationally ice models. It is worth mentioning that on the CWS, due to its inherent crystallinity, on overwhelming and in practice unfeasible. Here, we opted to simulate the reactions on one surface reaction site is available. However, this is not the case for the ASW model, one reaction site, which is the one present in the cavity of the ASW model. This choice is which different adsorption sites can be identified, which are potential reaction sites. Ne because interstellar ices are materials that seem to present some degree of porosity [1,77]. ertheless, computing By proceeding all the in this way, reactions on our simulations allwith align these surface a more sitesadsorption realistic is computationally scenario ove whelming and in practice unfeasible. Here, we opted to simulate the reactions on o reaction site, which is the one present in the cavity of the ASW model. This choice cause interstellar ices are materials that seem to present some degree of porosity [1,7 proceeding in this way, our simulations align with a more realistic adsorption sce Int. J. Mol. Sci. 2023, 24, 16824 from an interstellar ice perspective, while significantly reducing the potential 8 of 25 for d adsorption sites (and the computational cost). As mentioned in Section 4.2, for the calculations with these periodic system from adoptedan interstellar ice perspective, the cost-effective while significantly DFT//HF-3c approach reducing the potential to simulate for diverse in whic the reactions, adsorption sites (and the computational cost). DFT functional is that chosen in the benchmarking study completed in the gas phas As mentioned in Section 4.2, for the calculations with these periodic systems, we Appendix adopted theAcost-effective for more details on theapproach DFT//HF-3c choice of tothe functional). simulate For all the reactions, the reactions, in which the th step DFT was the adsorption functional ofinthe is that chosen thereactants by placing benchmarking them manually study completed in the gas on the(see phase surface. Appendix A for more details on the choice of the functional). For all the reactions, the first step was 2.2.1. the adsorption Formamide of the reactants by placing them manually on the surface. Formation 2.2.1.The formation Formamide pathway towards formamide (see reaction 1) was simulated o Formation CWSThe and ASW models, formation pathwayadopting the M06-2X-D3//HF-3c towards formamide (see Reaction (1))methodology. was simulated onFigure the 5 dis CWS and ASW models, adopting the M06-2X-D3//HF-3c methodology. Figure 5 displays the optimized geometries of the stationary points of the reaction and the ZPE-corr the optimized PES simulatedgeometries of the on the two stationary surface points of the reaction and the ZPE-corrected models. PES simulated on the two surface models. Figure 5. M06-2X-D3//HF-3c ZPE-corrected PES and the corresponding optimized stationary points, Figure 5. M06-2X-D3//HF-3c ZPE-corrected PES and the corresponding optimized stationary p alongside the relative energies, for the formamide formation on the CWS (A) and ASW (B) models. alongside the relative energies, for the formamide formation on the CWS (A) and ASW (B) m The unit cell is highlighted in black in the CWS model. Distances are in Å and relative potential The unit cell is −highlighted in black in the CWS model. Distances are in Å and relative po energies in kJ mol 1 . Color coding: white, H atoms; grey, C atoms; blue, N atoms; and red, O atoms. energies in kJ mol . Color coding: white, H atoms; grey, C atoms; blue, N atoms; and red, O −1 The thermodynamics of the reaction remain consistent with those observed in th phase. The product, which includes HCONH2 and H adsorbed on the surface, repre Int. J. Mol. Sci. 2023, 24, 16824 9 of 25 The thermodynamics of the reaction remain consistent with those observed in the gas phase. The product, which includes HCONH2 and H adsorbed on the surface, represents the most stable stationary point of the energy profiles. However, it is noteworthy that the difference in relative energies between the product and the intermediate is more pronounced in both surface models (18.2 kJ mol−1 for CSW and 8.3 kJ mol−1 for ASW) than in the gas-phase results (3.7 kJ mol−1 ). This is because they are stabilized by favorable interactions with the surface, such as H-bonds and dispersion forces, which are absent in the gas-phase scenario. Regarding the first step of the reaction (the coupling between H2 CO and NH2 ), it presents very high energy barriers (63.7 and 51.7 kJ mol−1 on CWS and ASW, respectively), which contrasts to that in the gas phase (22.6 kJ mol−1 ), probably due to the need to break the reactant/surface interactions to reach the transition structure (which is not necessary in the gas phase). The difference between the energy barriers when on CWS or ASW models is probably due to the high rigidity of the CSW model, which inhibits stabilizing surface reconstructions that are possible in the softer ASW model. The second step (dissociation of the intermediate), in contrast to the first step, presents energy barriers of 35.6 and 38.3 kJ mol−1 on CWS and ASW, which are significantly lower than that in the gas phase (45.4 kJ mol−1 ). 2.2.2. Acetaldehyde Formation The acetaldehyde formation pathway (Reactions (2) and (3)) was simulated on the water ice models through employing the BHLYP-D3(BJ)//HF-3c methodology. In the benchmark study, the M06-2X-D3 functional provided the most accurate description of these reactions, but we encountered dramatic energy convergence problems in the single- point calculations within the periodic approach. Consequently, we moved on the BHLYP- D3(BJ) method, as it was the second most accurate method describing the reactions. The geometries of the stationary points of the reactions and the ZPE-corrected PES simulated on the CWS and ASW models are depicted in Figure 6. Regarding the first reaction, the products formed on both surfaces (CH3 CHOH + OH) constitute the most stable stationary point of the profile, as also found in the gas-phase scenario. In contrast, the identified energy barriers for this reaction are 5.6 kJ mol−1 on the CWS and barrierless on the ASW (in this latter case when ZPE-corrections are considered), which are significantly lower than that in the gas phase (18.1 kJ mol−1 ). As far as the second reaction is concerned, the results obtained on the CWS and ASW follow the same trend as in the gas phase, that is, a barrierless O addition to CH3 CHOH forming the CH3 CHOOH intermediate, which is the most stable structure of the energy profile, meaning that the final acetaldehyde product is thermodynamically disfavored because it is less stable than the intermediate. Moreover, the dissociation step of the intermediate to form acetaldehyde presents energy barriers of 341.8 and 217.8 kJ mol−1 on CWS and ASW, respectively, which are notably higher than that computed in the gas phase (72.8 kJ mol−1 ). The energy barrier difference between the CWS and the ASW can be attributed to the highest rigidity of the CWS that hampers the stabilizing surface reconstruction effects that are possible to take place in the less-rigid ASW system. 2.2.3. Methyl Formate Formation The formation of methyl formate adopting the Reaction (4) was simulated on the CWS and ASW models at the B3LYP-D3(BJ)//HF-3c theory level. Figure 7 displays the geometries of the stationary points of the reaction and the ZPE-corrected PES simulated on the two water ice surface models. Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 10 of 26 Int. J. Mol. Sci. 2023, 24, 16824 10 of 25 Figure 6. BHLYP-D3(BJ)//HF-3c ZPE-corrected PES and the corresponding optimized stationary Figure 6. BHLYP-D3(BJ)//HF-3c ZPE-corrected PES and the corresponding optimized stationary points, alongside the relative energies, for the acetaldehyde formation on the CWS (A) and ASW (B) points, alongside the relative energies, for the acetaldehyde formation on the CWS (A) and ASW (B) models. The unit cell is highlighted in black in the CWS model. Distances are in Å and relative models. The unit cell is highlighted in black in the CWS model. Distances are in Å and relative potential energies in kJ mol−1 .−1Color coding: white, H atoms; grey, C atoms; and red, O atoms. potential energies in kJ mol . Color coding: white, H atoms; grey, C atoms; and red, O atoms. Regarding the first reaction, the products formed on both surfaces (CH3CHOH + OH) constitute the most stable stationary point of the profile, as also found in the gas-phase scenario. In contrast, the identified energy barriers for this reaction are 5.6 kJ mol−1 on the 2.2.3. Methyl Formate Formation The formation of methyl formate adopting the reaction 4 was simulated on the CWS and ASW models at the B3LYP-D3(BJ)//HF-3c theory level. Figure 7 displays the geome- Int. J. Mol. Sci. 2023, 24, 16824 11 of 25 tries of the stationary points of the reaction and the ZPE-corrected PES simulated on the two water ice surface models. 7. B3LYP-D3(BJ)//HF-3c Figure 7. Figure B3LYP-D3(BJ)//HF-3c ZPE-corrected ZPE-correctedPES, PES,and andthethe corresponding corresponding optimized stationary optimized stationary points, alongside points, alongsidethe therelative relativeenergies, for for energies, the the methyl formate methyl formation formate on the on formation CWSthe(A) and(A) CWS ASW (B)ASW and (B) models. models. TheThe unitunit cellcell is highlighted is highlighted in black in black in CWS in the the CWS model. model. Distances Distances are are in Å in Å relative and and relative potential energiesin potential energies inkJ mol−−11 .. Color kJmol Color coding: coding:white, white,HHatoms; atoms; grey, grey, C atoms; C atoms; andand red,red, O atoms. O atoms. On the water ice models, the O addition to the CH2 OCH3 radical to form the H2 COOCH3 intermediate is also barrierless, as occurred in the gas phase, due to being a spin–spin- coupling driven reaction. Dissociation of this intermediate gives rise to the final reaction. Methyl formate formation pathway presents an energy barrier of 32.5 kJ mol−1 on the CWS model, which is significantly lower than in the gas phase (49.0 kJ mol−1 ). However, on the ASW, this step presents an energy barrier as high as 97.9 kJ mol−1 (higher than in the gas phase). The elevated energy barrier in the ASW model is attributed to a substantial internal structural reorganization of the surface during the dissociation process, giving rise to this elevated energy barrier, while on the CWS model, the dissociation step occurs without significant structural reorganization. Regarding the stability of the final products, on both CWS and the ASW, they present a similar stability as the intermediate; on the CWS the products are just 3.3 kJ mol−1 higher in energy, while on the ASW model they are just 2.0 kJ mol−1 lower in energy. This trend is analogous to the gas-phase reaction, where both the intermediate and the products display similar stability, the intermediate being the most stable structure of the profile and 6.7 kJ mol−1 more stable than the products. Int. J. Mol. Sci. 2023, 24, 16824 12 of 25 2.2.4. Formic Acid Formation The formation of formic acid, through the Reaction (5), was simulated on the CWS and ASW Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW models using the B3LYP-D3(BJ)//HF-3c approach. The ZPE-corrected PES and 13 of 26 the geometries of the stationary points involved in this formation route on the two water ice models are depicted in Figure 8. Figure 8. B3LYP-D3(BJ)//HF-3c ZPE-corrected PES and the corresponding optimized stationary Figure alongside points, 8. B3LYP-D3(BJ)//HF-3c ZPE-corrected the relative energies, PES and for the formic acid the corresponding formation optimized on the CWS (A) andstationary ASW (B) models. The unit cell is highlighted in black in the CWS model. Distances are in(A) points, alongside the relative energies, for the formic acid formation on the CWS Å and and ASW (B) relative models. The unit cell is highlighted − 1 in black in the CWS model. Distances are in Å and relative potential energies in kJ mol . Color coding: white, H atoms; grey, C atoms; and red, O atoms. potential energies in kJ mol−1. Color coding: white, H atoms; grey, C atoms; and red, O atoms. The results obtained on CWS and ASW indicate that the presence of the water ices Thethe reverses results obtained energetic onwith trend CWSrespect and ASW indicate to the that the gas-phase presence scenario. of the That water is, on the ices icy reverses the energetic trend with respect to the gas-phase scenario. That is, on the surfaces, the final products (HCOOH and H on the surfaces) are the most stable stationary icy sur- faces, of points thethe final products path. (HCOOH This was and Hinonthe also observed the surfaces)ofare formation the most on formamide stable stationary the water ices points of the path. This was also observed in the formation of formamide on the water ices (see above) and can be attributed to the emergence of favorable interactions between the newly formed products and the surface, such as H-bond and dispersion, which are miss- ing in gas-phase conditions. The coupling step presents high energy barriers on the water ice surface models (59.5 Int. J. Mol. Sci. 2023, 24, 16824 13 of 25 (see above) and can be attributed to the emergence of favorable interactions between the newly formed products and the surface, such as H-bond and dispersion, which are missing in gas-phase conditions. The coupling step presents high energy barriers on the water ice surface models (59.5 and 83.6 kJ mol−1 on CWS and ASW, respectively), which are significantly higher energy barriers than those in the gas phase (33.1 kJ mol−1 ), due to the need to break the reactants/surface interactions. In contrast, the dissociation step presents energy barriers of 23.7 and 22.1 kJ mol−1 on CWS and ASW, respectively, which are lower than those obtained in the gas phase (59.2 kJ mol−1 ). 3. Discussion In this work, we modeled the formation pathways of four astrochemically relevant iCOMs: formamide, acetaldehyde, methyl formate, and formic acid. As formation routes, we took feasible gas-phase neutral–neutral reactions, which, additionally, were also sim- ulated on models of CWS and ASW ices. Table 1 summarizes the energetics of the simu- lated reactions. Table 1. Computed ZPE-corrected energetics (in kJ mol−1 ) of the reactions simulated under different conditions (gas phase, and on the CWS and the ASW models). ∆E‡ stands for the intrinsic energy barriers and ∆EX−Y for the relative energy between X and Y. Reaction Conditions Energy Terms ∆E1 ‡ ∆E2 ‡ ∆EI−R ∆EP−R Gas phase 22.6 45.6 −34.4 −38.1 Reaction (1) CWS 63.7 35.6 −41.6 −59.8 ASW 51.7 38.4 −36.7 −45.0 ∆E‡ ∆EP1−R1 Gas phase 18.1 −95.0 Reaction (2) CWS 5.6 −69.5 ASW Barrierless −89.8 ∆E1 ‡ ∆E2 ‡ ∆EI−R2 ∆EP2−R2 Gas phase Barrierless 72.8 −364.9 −312.0 Reaction (3) CWS Barrierless 341.8 −331.8 −299.2 ASW Barrierless 271.8 −319.4 −259.3 ∆E1 ‡ ∆E2 ‡ ∆EI−R ∆EP−R Gas phase Barrierless 49.0 −381.6 −374.9 Reaction (4) CWS Barrierless 32.6 −365.9 −362.5 ASW Barrierless 98.0 −379.2 −381.1 ∆E1 ‡ ∆E2 ‡ ∆EI−R ∆EP−R Gas phase 33.1 59.2 −79.6 −63.1 Reaction (5) CWS 59.5 23.7 −61.6 −71.9 ASW 83.6 22.1 −43.1 −57.2 All these processes involve multi-step mechanisms. In the gas-phase scenario, the reaction intermediates are usually the most stable stationary points, and when they are not, are very similar in energy as the final products. This could lead us to think that the processes are thermodynamically disfavored, but (as mentioned in the Introduction), in the gas phase, whether a reaction is favorable or not is not only regulated by their energetic features, but overall kinetics need to be accounted for, normally by means of calculations based on the RRKM theory. In this sense, available studies indicate that these paths are actually favorable and hence they have been proposed to be major sources for these iCOMs. A particular feature enabling the occurrence of these reactions is that the energy of the transition states is lower than the asymptotes, meaning that the system has enough internal energy overcoming the energy barriers. On the other hand, when these reactions are simulated on the periodic water ice models, in some cases, the results indicate that the Int. J. Mol. Sci. 2023, 24, 16824 14 of 25 presence of the surface alters the energetic trends compared to the gas-phase conditions. That is, in the formation of formamide and formic acid, the reaction products are identified as the most stable structures, more than the respective intermediates. In the formation of methyl formate, the stabilities between the products and the intermediate are very similar, while the formation of acetaldehyde is the unique case in which the intermediate is still the most favorable structure along the path. Modification of the thermodynamics of the reactions on the surface with respect to the gas-phase scenario is due to the interactions between the structures and the surface, which cannot take place in the gas phase. These interactions (mostly based on H-bonding and dispersion) have been found to be indeed important in driving reactivity of similar organic molecules, which are used as conformers in co-crystals in relevant processes for pharmaceutical industry [78–81]. Although all the reactions considered in this work are neutral–neutral reactions, to facilitate the discussion, we can categorize the steps studied as (i) H abstraction, which is the first reaction of the formation mechanism of acetaldehyde (Reaction (2)); (ii) addition of atomic O to a radical, the first step of the formation mechanism of methyl formate (Reaction (4)) and also the first step of the second reaction of the formation mechanism of acetaldehyde (Reaction (3)); (iii) coupling of formaldehyde with a radical, in particular with NH2 (first step of formamide formation, Reaction (1)) and with OH (first step of formic acid formation, Reaction (5); and (iv) dissociation of intermediates, which occurs in the second step of formamide, methyl formate, and formic acid formations (Reactions (1), (4), and (5), respectively) and also in the second step of the second reaction of acetaldehyde formation (Reaction (3)). When considering the reaction belonging to the category (i)—H abstraction—in the gas phase, it exhibits a high intrinsic energy barrier (18.1 kJ mol−1 ). As mentioned above, despite this barrier, the TS is 1.3 kJ mol−1 below the asymptote (see Figure 2) and hence the reaction is assumed to be favorable under these conditions. When this process is simulated on the periodic water ice models, a transition state is identified on both models. On the CWS, the ZPE-corrected barrier is 5.6 kJ mol−1 , while on the ASW, after adding the ZPE-corrections, the barrier is submerged below the reactants, and accordingly it can be considered barrierless. Irrespective of the model, it is clear that the presence of the water ice surface decreases the energy barriers, and this way exerting a catalytic effect. This is probably due to the interactions of the CH2 moiety of ethanol and the OH radical with the surface which polarize the C-H and O-H chemical bonds, facilitating the H abstraction. Therefore, considering that interstellar water ices are amorphous, we can state that these H abstraction processes are likely to occur in the ISM as a grain-surface reaction. When considering the processes belonging to the category (ii)—O addition to radicals –either in the gas phase or on the water ice surfaces, they are barrierless. This is because the process is driven by a spin–spin coupling, which under gas-phase conditions are definitely barrierless, but it could not be the case on water ice surfaces. Indeed, other radical–radical processes (also driven by spin–spin couplings) have shown that these processes can present energy barriers, the emergence of which are primarily attributed to the interactions between the radicals and the surface (mainly H-bonds and dispersion), which need to disrupt so that the reactions proceed. In a study by some of us [29], it was established that the interactions with the surface play a crucial role in the magnitude of these energy barriers. Specifically, strong radical–surface interactions (like H-bonds) are expected to result in relatively high energy barriers because it is necessary to break them for the progress of the reaction, which is translated to an increased energy barrier. In contrast, for weak radical-surface interactions (like the dispersive forces), low energy barriers or even barrierless reactions are expected, because reaching the transition state structure involves essentially the rotation and/or translation of the reactive species, which are low energetic motions. In our specific cases, both the CH3 CHOH and CH3 OCH2 radicals are firmly attached to the water ice surfaces due to forming H-bonds. However, atomic O is weakly bounded to the surface (it indeed presents a very low binding energy, 13 kJ mol−1 [82], enabling an efficient diffusion on icy water surfaces [83,84]). In our simulations, we observed that only a slight rotation of the Int. J. Mol. Sci. 2023, 24, 16824 15 of 25 radical to become well oriented with the O atom was required to proceed with the coupling, which resulted in barrierless reactions. Therefore, these processes can also be considered as feasible in the ISM. When considering the processes falling into category (iii)—formaldehyde–radical couplings—we observe an increase of the energy barriers when they occur on the water ice surfaces when compared with the gas-phase ones. This can be attributed to the favorable interactions between the reactants and the surface leading to an extra-stabilization of the reactants, which is not possible in the gas phase. However, we did not find a consistent trend between the energy barriers on the CWS and the ASW models. That is, in the formamide case the energy barrier is larger on CWS than on ASW, while in the formic acid case it is the opposite. As mentioned above, to enable the system to overcome the activation barriers, it is required to break the interactions between reactants and the surface, and thus stronger interactions with the surface are predicted to result in higher activation barriers. It is worth mentioning that in the works of Perrero et al. [85], in which the binding energies (BEs) of several species using the same periodic ice models as here were computed, it was found that a distribution of BEs could be determined when using the ASW model due to the presence of different binding sites, which are also different reaction sites. Therefore, depending on the initial positions of the reactants and their interactions with the surface, the energy barrier can either increase or decrease when transitioning from the CWS model to the ASW model. Irrespective of that, these formaldehyde–radical coupling processes present substantial energy barriers and accordingly they could hardly happen in interstellar conditions. When considering processes of category (iv)—intermediate dissociation—they present substantial energy barriers when occurring both in the gas-phase conditions (considering intrinsic energy barriers) and on the water ice models. As mentioned above, neutral– neutral reactions usually lead to the formation of a radical intermediate, which under strict gas-phase conditions cannot be stabilized by dispersing its internal energy, and, thus, the reaction can dissociate back or evolve to form new species (in exothermic reactions). In contrast, when these reactions occur on water ice surfaces, the water ices can act as third bodies (as demonstrated by some of us in recent works [19,86,87]), in which the intermediate is stabilized by transferring the energy released by the reaction to the ice. Consequently, the intermediate becomes stable on the surface and does not contain the same excess of internal energy as occurs in the gas phase, and accordingly substantial dissociation energy barriers emerge, which are hard to surmount in ISM conditions. However, when the intermediate dissociation is associated with the release of a H atom (as is the case of formamide, methyl formate, and formic acid formation), tunneling effects can become relevant, especially at very low temperatures [88,89]. Tunneling is more likely in symmetric processes than in non- symmetric processes, that is, when the initial and final states present similar stabilities. This is the case of methyl formate on the water ices, in which the intermediate and its dissociated product are almost isoenergetic. To check upon the possibility of tunneling, we conducted a simplistic (but useful) assessment by calculating the crossover temperature (TX ), namely, that temperature below which tunneling dominates, by adopting the Fermann Auerbach formulation [90], which accounts directly for the energy barrier height and indirectly for the energy barrier width with the transition frequency. The results show, for this process, TX of 217 and 399 K on CWS and ASW, respectively, and accordingly tunneling can dominate for the interstellar conditions, making this step possibly favorable. Bearing in mind that the first step of the methyl formate formation is the O addition to CH2 OCH3 (barrierless), we can state that its formation through Reaction (4) is plausible in the ISM. In the formation of formamide and formic acid, the intermediate dissociation is a non-symmetric process, in which the final product is more stable than the intermediate, which is also a favorable situation for the dominance of tunneling (but less than for symmetric barriers). To assess the tunneling effects on these two processes, we also calculated the TX , which are 233 and 233 K, and 223 and 257 K on CWS and ASW in the formamide and formic acid formation, respectively. Thus, the intermediate dissociation by releasing an H atom in the final step of Int. J. Mol. Sci. 2023, 24, 16824 16 of 25 the formamide and formic acid formation seems to be feasible under tunneling regimes in the ISM. Nevertheless, it is worth remembering that, at variance with the methyl formate formation, the initial steps in the formation routes of these two compounds have been found to be unlikely in the ISM due to presenting high energy barriers, so that the formation of these iCOMs is highly hampered. Despite this, it is worth mentioning that the presence of the water ice surfaces modifies the relative energies between the intermediates and the products, always in favor of the stabilization of the products. Thus, the icy surfaces allow the obtaining of energy barrier situations (namely, symmetric, and exoenergetic non-symmetric) favorable for the occurrence of tunneling. Obviously, a more rigorous treatment of the tunneling effects is needed to definitely determine their actual role in these steps, such as applying semi-classical approaches to the transition state and RRKM theories [91] or the instanton theory [92]. In summary, our results predict that gas-phase processes can hardly occur on water ice surfaces. The only reaction which can potentially occur is the formation of methyl formate, in which tunneling effects need to be advocated. Thus, a main message of this work is that there is a clear frontier between gas-phase and grain-surface chemistry, that is, gas-phase processes can be feasible as specifically gas-phase reactions but not as grain- surface reactions, and vice versa. Consequently, results emerging from one scenario cannot be directly translated to the other. Therefore, establishing a general trend for the lack of analogy between gas-phase and grain-surface mechanisms is challenging, as the feasibility of the reactions depends on the nature of the reaction steps, the reactive species involved, and the surface model, thus requiring a case-by-case investigation. 4. Materials and Methods 4.1. Surface Modelling A periodic approach was adopted to model the surfaces mimicking the interstellar water ice mantles. Two different models were considered: one crystalline and another one amorphous. The crystalline ice surface model was based on the P-ice structure. P-ice represents the proton ordered counterpart to the hexagonal (proton-disordered) water ice, which has been shown to replicate well the physico-chemical properties of crystalline water ice [93]. The 2D periodic surface model was generated by cutting the P-ice 3D periodic bulk structure perpendicular to the [010] direction. As a result, the (010) slab surface model was obtained, which is one of the most stable planes among crystalline P-ice surfaces, presenting a null dipole component perpendicular to the surface [94,95]. The resulting crystalline water surface model (hereafter referred to as CWS) contains 12 atomic layers, with a total thickness of 10.90 Å, and includes 24 water molecules (72 atoms) per unit cell. However, the CWS model is not representative of the structural features of actual interstellar water ices, as they are predominantly found in an amorphous state [1,96]. To better depict a more realistic interstellar water ice surface, an amorphous model designed by Ferrero et al. [97] was also employed. The amorphous surface was constructed by assembling small clusters generated by Shimonishi et al. [98]. Like the CWS model, this amorphous solid water model (hereafter referred to as ASW) is based on a periodic ap- proach but in order to account for the amorphousness of the system, a larger unit cell was used (20.35 Å of length), containing 60 water molecules (180 atoms) per unit cell. Figure 9 displays the top and lateral view of CWS and the lateral view of ASW structural models. Int. J. Mol. Sci. 2023, 24, x FOR PEER REVIEW 18 of 26 Int. J. Mol. Sci. 2023, 24, 16824 17 of 25 Figure 9. Figure Top view 9. Top view (A) (A) and and lateral lateral view view(along (alongthe thebblattice latticevector) vector)(B) (B)ofofthe thecrystalline crystallinewater water surface surface model. Lateral model. Lateral view view (along (along the the aa lattice lattice vector) vector)ofofthe theamorphous amorphoussolid solidwater watersurface surface(C). Unit (C). cellcell Unit highlighted in black. Color highlighted Color coding: coding: white, white,H Hatoms; atoms;andandred, red,OOatoms. atoms. 4.2. Computatonal 4.2. Computatonal DetailsDetails Simulations were Simulations were carried carried out out employing employingthe theperiodic periodicab abinitio initioCRYSTAL17 CRYSTAL17code code[99]. [99]. This software implements both the Hartree–Fock and This software implements both the Hartree–Fock and the Kohn–Sham self-consistentthe Kohn–Sham self-consistent field field methods to methods to solve solve thetheelectronic electronicSchrödinger Schrödinger equation. equation. CRYSTAL17 CRYSTAL17 usesuseslocalized Gaussian localized Gauss- type orbitals (GTOs) as basis sets, which are functions centered ian type orbitals (GTOs) as basis sets, which are functions centered to the atoms that decay to the atoms that decay with the distance with respect to each nucleus, allowing accordingly the simulation of both with the distance with respect to each nucleus, allowing accordingly the simulation of periodic (crystals, surfaces, and polymers) and non-periodic (molecules). This is at variance both periodic (crystals, surfaces, and polymers) and non-periodic (molecules). This is at with basis sets based on plane wave (PWs) functions, which are periodic functions that fill variance with basis sets based on plane wave (PWs) functions, which are periodic func- all the 3D space. Consequently, the definition of surfaces in CRYSTAL17 does not require tions that fill all the 3D space. Consequently, the definition of surfaces in CRYSTAL17 does the 3D replica commonly used in PW-based codes, and thus they are represented as true not slabrequire models,the with 3Da replica complete commonly empty space used in PW-based above and belowcodes, and thus they are repre- the slab. sented The characterization of the potential energy surface (PES) of thebelow as true slab models, with a complete empty space above and reactionsthe required slab. The characterization determines the structures and of thethepotential energeticsenergy of theirsurface stationary(PES) of the points, reactions including required reactants, determines the structures and the energetics of their stationary intermediates, products, and transitions states. A cost-effective approach was employed for points, including reactants, intermediates, this purpose, referredproducts, to asand transitionsInstates. DFT//HF-3c. A cost-effective this composite procedure, approach was employed the geometries of the for this purpose, stationary points were referred to as DFT//HF-3c. optimized at HF-3c level, In and, this composite subsequently, procedure, an energy therefinement geometries of thewas step stationary carried points were optimized out through single-pointatenergy HF-3c calculations level, and, subsequently, with a densityan energy re- functional finement theory (DFT) step method was carried on the outoptimized through single-point HF-3c geometries.energy HF-3c calculations with a density is a semi-empirical method within functional theory the Hartree–Fock (DFT) method on framework, the optimized using a minimal HF-3c Gaussian geometries. HF-3c atomic orbital is a semi-em- pirical method within the Hartree–Fock framework, using a minimal Gaussian atomictoor- (Gaussian-AO) basis set known as MINIX. It includes three a posteriori corrections alleviate bital the deficits arising (Gaussian-AO) basis setfrom the approximations known introduced. as MINIX. It includes These three correctionscorrections a posteriori account foralleviate to (i) short-ranged the deficits deficiencies arising from due to thethe use of a small or approximations minimal basis introduced. These set;corrections (ii) basis setac- superposition count error (BSSE)deficiencies for (i) short-ranged implementing duethe to geometrical the use of a counterpoise small or minimal correction basis gCP); set; (ii) and (iii) dispersive interactions [100–104]. basis set superposition error (BSSE) implementing the geometrical counterpoise correc- For the tion gCP); andrefinement of the interactions (iii) dispersive energetics, different [100–104]. DFT functionals were considered: the pure gradient generalized approximation (GGA) PBE; the non-local hybrids B3LYP and For the refinement of the energetics, different DFT functionals were considered: the BHLYP, which include a 20% and a 50% of exact exchange in its definition, respectively; and pure gradient generalized approximation (GGA) PBE; the non-local hybrids B3LYP and the meta-hybrid M06-2X, which incorporates a 45% of exact exchange [105–107]. The basic BHLYP, which include a 20% and a 50% of exact exchange in its definition, respectively; set used for these DFT calculations was the 6-311+G(d,p) one. Furthermore, to account and the meta-hybrid M06-2X, which incorporates a 45% of exact exchange [105–107]. The for dispersive interactions, the corresponding Grimme’s D3 correction within the Becke– basic Johnsonset used scheme for(i.e., thesetheDFT D3(BJ)calculations correction)was wasthe 6-311+G(d,p) included one. Furthermore, when available; to ac- if not available, count the D3 correction was used [108,109]. For both CWS and ASW, the Hamiltonian matrix wasthe for dispersive interactions, the corresponding Grimme’s D3 correction within Becke–Johnson diagonalized inscheme 4k points (i.e., thefirst in the D3(BJ) correction) Brillouin was included when zone, corresponding available; to a shrinking if not factor available, of 2. Values theofD3 correction 7, 7, 7, 7, 14 for wastheused [108,109]. tolerances For both controlling CWS and Coulomb ASW, and the Hamiltonian the exchange series matrix was diagonalized were adopted in 4k points for all calculations. in choice The the firstofBrillouin zone, corresponding the numerical values controlling to a these shrink- ing factor of 2. Values of 7, 7, 7, 7, 14 for the tolerances controlling computational parameters is fully justified in Ferrero et al. [97]. For geometry optimization, Coulomb and the ex- change series were the convergence adopted criterium was − 7 forsetalltocalculations. 10 Hartree,The choicethat meaning of the the numerical self-consistent values con- fields trolling these computational (SCF) procedure was stopped parameters when the energy is fully justifiedbetween difference in Ferrero et al. cycles was[97]. For geome- smaller than try thisoptimization, value. the convergence criterium was set to 10−7 Hartree, meaning that the self- consistent fields (SCF) procedure was stopped when the energy difference between cycles was smaller than this value. Int. J. Mol. Sci. 2023, 24, 16824 18 of 25 To identify the DFT functional that described more accurately each chemical reac- tion, a gas-phase benchmarking study was conducted using the molecular GAUSSIAN16 program [110]. This study involved optimizing geometries using the four different DFT functionals (PBE, B3LYP, BHLYP, M06-2X) with a 6-311+G(d,p) basis set and the ener- gies refined by performing single-point energy calculations at the same DFT methods but employing a 6-311++G(2 df,2 pd) basis set. To assess the accuracy of the DFT meth- ods, the results were confronted with those obtained at the CCSD(T) level (considered the “gold-standard” method of quantum chemistry), which were obtained by performing single-point energy calculations on the DFT-optimized geometries using the Dunning aug-cc-pVTZ basis set [111] (hereafter referred to as CCSD(T)//DFT methodology). By comparing the energy results computed at full DFT level with the CCSD(T)//DFT ones, the most suitable DFT functional for each reaction was selected (see Table A1). Moreover, since periodic simulations were based on a DFT//HF-3c scheme, we performed an addi- tional benchmarking for the formamide gas-phase formation comparing this methodology with the CCSD(T)//DFT and CCSD(T)//HF-3c ones. The obtained results (shown in Tables A1 and A2) show very good agreement between these methodologies, pointing out the reliability of the DFT//HF-3c scheme for the periodic systems. To prevent surface deformations, periodic geometry optimizations were carried out by relaxing the internal atomic positions but by keeping fixed to the lattice parameters. Additionally, the spin-unrestricted formalism was adopted for open-shell systems. The na- ture of the stationary points of the reactions was characterized by calculating the harmonic frequencies at the HF-3c level; the outcome being minima for reactants, intermediates, and products, and first-order saddle points for transitions states. In CRYSTAL17, har- monic frequency calculation is conducted numerically at the Γ point by diagonalizing the mass-weighted Hessian matrix of the second-order energy derivatives with respect to atom displacements (i.e., central difference formula), in this case ±0.003 Å from the minimum along each Cartesian coordinate [112,113]. Here, in frequency calculations, the convergence criterium was set to 10−10 Hartree. With the harmonic frequency calculations, the vibra- tional zero-point energy corrections (ZPE) were obtained and applied to the corresponding potential-energy values. 5. Conclusions In this study, we investigated the means of quantum-chemical simulations the forma- tion pathways of four astrochemically relevant iCOMs: formamide, acetaldehyde, methyl formate, and formic acid. We explored these pathways by adopting neutral–neutral gas- phase synthetic routes (considered as major production sources of these molecules in the gas phase) and exporting them on water ice surface models. Our approach involved char- acterizing the potential energy surfaces for these mechanisms in the gas phase and on the periodic water ice models (crystalline and amorphous ones). We adopted a composite procedure, the DFT//HF-3c approach, in which we optimized geometries at the HF-3c level and then performed single-point energy calculations at the DFT level on the optimized HF-3c geometries. The DFT methods were benchmarked against CCSD(T) calculations in the gas-phase reactions. The results obtained indicate that on the water ice models, while some reactions exhibit lower intrinsic energy barriers compared to the gas-phase results, or even barrierless situations, each mechanism investigated included at least one high energy barrier step, which at low temperatures are unlikely to be overcome, thus hampering the formation of these molecules on grain surfaces. The presence of such high energy barriers is due to the third body effect of the icy surfaces, which by dissipating the internal energies of the intermediates stabilize them. The only exception in which the iCOM can be formed adopting a gas-phase neutral–neutral mechanism on the water ice surfaces is through methyl formate, in which the energetic step can be dominated by tunneling effects. An interesting outcome emerging from this work is that the presence of the surface can modify the overall energetics of the reactions compared with the gas-phase processes. Int. J. Mol. Sci. 2023, 24, 16824 19 of 25 While in this latter scenario, reaction intermediates are usually found to be more stable than the final products, on the water ice surface models this is not the case, in which the final products are in general more stable than the intermediates, thanks to the favorable interactions with the surface. In conclusion, this study highlights the fundamental differences between gas-phase and grain-surface chemistry. Formation routes cannot be translated between these sce- narios due to the modifying influence of the surface in the energy trends of the reactions. Accordingly, establishing a generalization between gas-phase and grain-surface chemistry is challenging, and, therefore, we advocate that a case-by-case study is required. Author Contributions: Conceptualization, A.R.; methodology, B.M.-B.; software, B.M.-B.; validation, B.M.-B. and A.R.; formal analysis, B.M.-B.; investigation, B.M.-B. and A.R.; resources, A.R.; data curation, B.M.-B. and A.R.; writing—original draft preparation, B.M.-B.; writing—review and editing, A.R.; visualization, B.M.-B. and A.R.; supervision, A.R.; project administration, A.R.; and funding acquisition, A.R. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement No. 865657 for the project “Quantum Chemistry on Interstellar Grains” (QUANTUMGRAIN), and by MICINN for the project PID2021-126427NB-I00. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data are contained within the article and in Appendix A. Acknowledgments: We thankfully acknowledge the computer resources and assistance provided by the Barcelona Supercomputing Center (BSC-RES) and CSUC. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. Appendix A. Benchmark Study To select the DFT function that most accurately describes each reaction to perform an energy refinement in the periodic approach, we initially conducted a benchmark study in the gas-phase, specifically in the absence of water molecules. Geometry optimizations of the stationary points of the reactions were performed using four different DFT functionals. To account for dispersion forces, we included the Grimme’s D3(BJ) dispersion term when available; otherwise, the D3 correction was used. In particular, the DFT functionals used were PBE-D3(BJ), B3LYP-D3(BJ), BHLYP-D3(BJ) and M06-2X-D3, in conjunction with the 6-311+G(d,p) basis set for geometry optimizations and the 6-311++G(2 df,2 pd) basis set for energy refinement. Additionally, we refined the energies by performing single-point energy calculations at CCSD(T) level of theory using the Dunning aug-cc-pVTZ basis set. The DFT results were confronted with the refinement at CCSD(T) level of theory. Table A1 reports the relative energies with respect to the reactants of the stationary points of each reaction using the four different DFT methods, along with the refinement at the CCSD(T) level. The first step of the second reaction for the formation of acetaldehyde and the first step for the formation of methyl formate are radical–radical couplings, and, therefore, they are barrierless. Significant differences are observed in the relative energies between the four different DFT functionals, although changes in geometries are almost insignificant, except for PBE-D3(BJ). PBE-D3(BJ) results indicates that molecule–radical reactions (first steps of formamide and acid formic reactions) are barrierless while all the other DFT functionals indicate the presence of a transition state. The benchmark study indicates that the best DFT for formamide and acetaldehyde reactions, is the M06-2X-D3 functional, while for methyl formate and formic acid reactions, it is the B3LYP-D3(BJ) functional. Due Int. J. Mol. Sci. 2023, 24, 16824 20 of 25 to SCF convergence problems in the periodic calculations for acetaldehyde reaction with M06-2X-D3 functional, we adopted the second more accurate DFT functional for periodic simulations, which is BHLYP-D3(BJ). To assess the accuracy of the DFT//HF-3c procedure applied in the periodic approach, a second benchmark study in the gas-phase was conducted for formamide reaction. In this study, geometries were optimized at HF-3c level of theory. Subsequently, single-point energy calculations at both DFT level and CCSD(T) level of theory were performed, using the same basis set employed in the first benchmark study. The DFT results were confronted with the refinement at CCSD(T) level of theory. Table A2 presents the relative energies with respect to the reactants for the stationary points of the reaction, employing the four different DFT functionals along with the refinement at the CCSD(T) level. The obtained results align with the findings of the previous benchmark study, identifying the M06-2X DFT functional as the most accurate for describing the formamide reaction. Table A1. Relative potential energies of the stationary points of the reactions under gas-phase condi- tions, for the DFT-optimized systems and for the CCSD(T)//DFT single-point energy calculations. Energy units are in kJ mol−1 . Relative errors of the full DFT values with respect to CCSD(T)//DFT for the energy barriers are also included (∆Ei ‡ (%)). PBE-D3(BJ) B3LYP-D3(BJ) BHLYP-D3(BJ) M06−2X-D3 Reaction Species PBE-D3(BJ) CCSD(T) B3LYP-D3(BJ) CCSD(T) BHLYP-D3(BJ) CCSD(T) M06−2X-D3 CCSD(T) H2 CO + NH2 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 H2 CO − NH2 – – −13.6 −10.9 −13.4 −11.1 −16.5 −11.0 TS1 – – −8.6 2.1 3.7 2.6 −1.6 3.3 H2 CONH2 −103.2 −62.6 −77.4 −63.2 −77.9 −63.7 −79.5 −63.7 HCONH2 TS2 −52.9 4.3 −14.9 5.4 4.3 5.1 −15.0 6.4 HCONH2 + H −75.5 −37.1 −42.1 −37.2 −31.7 −37.4 −57.7 −37.1 Error ∆E1 ‡ (%) – – 61.6 – 25.1 – 3.8 – Error ∆E2 ‡ (%) 24.8 – 8.9 – 19.6 – 7.9 – Average error (%) – – 35.3 – 22.4 – 5.9 – CH3 CH2 OH + OH 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 CH3 CH2 OH − OH −33.1 −29.1 −31.6 −29.4 −32.0 −29.5 −28.3 −26.8 TS1 −28.1 −21.7 −18.9 −8.4 6.7 2.0 −2.9 0.4 CH3 CHOH + H2 O −118.5 −92.6 −100.4 −92.4 −80.3 −92.2 −98.2 −92.5 CH3 CHOH + O 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 CH3 CHO CH3 CHOOH −457.3 −393.7 −412.1 −398.6 −356.9 −379.7 −394.7 −378.2 TS2 −348.7 −323.2 −315.8 −297.4 −261.1 −297.3 −296.6 −296.9 CH3 CHO + OH −328.9 −307.1 −313.9 −307.1 −279.4 −306.1 −308.4 −306.8 Error ∆E1 ‡ (%) 32.3 – 39.4 – 22.8 – 6.7 – Error ∆E2 ‡ (%) 55.3 – 4.9 – 16.3 – 20.8 – Average error (%) 43.8 – 22.1 – 19.5 – 13.8 – CH3 OCH2 + O 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 H2 COOCH3 −458.4 −378.6 −403.4 −393.5 −368.5 −394.2 −434.1 −394.5 HCOOCH3 TS −395.6 −328.7 −337.1 −327.3 −275.3 −325.3 −361.8 −325.7 HCOOCH3 + H −414.2 −365.3 −359.7 −365.4 −306.3 −363.8 −399.9 −365.0 Error ∆E‡ (%) 25.7 – 0.3 – 35.4 – 5.1 – H2 CO + OH 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 H2 CO − OH – – −23.1 −22.0 −23.9 −23.4 −23.2 −23.1 TS1 – – −5.3 7.8 12.9 7.9 −17.0 8.6 H2 COOH −154.8 −100.1 −109.1 −95.6 −116.7 −107.7 −126.1 −107.2 HCOOH TS2 −93.5 −38.5 −33.7 −20.6 −24.4 −36.5 −53.7 −36.0 HCOOH + H −113.2 −76.6 −56.5 −58.9 −55.7 −75.8 −92.0 −76.4 Error ∆E1 ‡ (%) – – 40.1 – 17.8 – 80.5 – Error ∆E2 ‡ (%) 0.6 – 0.6 – 29.6 – 1.6 – Average error (%) – – 20.4 – 23.7 – 41.1 – Int. J. Mol. Sci. 2023, 24, 16824 21 of 25 Table A2. Relative potential energies of the stationary points for the reaction of formamide formation under gas-phase conditions, calculated at DFT/HF-3c and CCSD(T)//HF-3c. Energy units are in kJ mol−1 . Relative errors with respect to CCSD(T)//HF-3c for the energy barriers are also included (∆Ei ‡ (%)). Species PBE-D3(BJ) B3LYP-D3(BJ) BHLYP-D3(BJ) M06−2X-D3 CCSD(T) H2 CO + NH2 0.0 0.0 0.0 0.0 0.0 H2 CO − NH2 −9.0 −5.6 −2.7 −4.9 −7.2 TS1 −19.1 4.4 22.8 16.5 16.2 H2 CONH2 −80.8 −61.2 −68.1 −67.8 −50.5 TS2 −46.3 −7.8 14.0 −8.2 6.8 HCONH2 + H −74.5 −41.3 −30.8 −57.2 −36.2 Error ∆E1 ‡ (%) 143.0 57.1 8.7 8.5 – Error ∆E2 ‡ (%) 39.9 6.8 43.4 4.0 – Average error (%) 91.4 32.0 26.0 6.3 – References 1. Boogert, A.A.; Gerakines, P.A.; Whittet, D.C. Observations of the icy universe. Ann. Rev. Astron. Astrophys. 2015, 53, 541–581. [CrossRef] 2. McClure, M.K.; Rocha, W.R.M.; Pontoppidan, K.M.; Crouzet, N.; Chu, L.E.U.; Dartois, E.; Lamberts, T.; Noble, J.A.; Pendleton, Y.J.; Perotti, G.; et al. An Ice Age JWST inventory of dense molecular cloud ices. Nat. Astron. 2023, 7, 431–443. [CrossRef] 3. McGuire, B.A. 2021 Census of Interstellar, Circumstellar, Extragalactic, Protoplanetary Disk, and Exoplanetary Molecules. Astrophys. J. Suppl. Ser. 2022, 259, 30. [CrossRef] 4. Herbst, E.; van Dishoeck, E.F. Complex Organic Interstellar Molecules. Ann. Rev. Astron. Astrophys. 2009, 47, 427–480. [CrossRef] 5. Herbst, E. The synthesis of large interstellar molecules. Int. Rev. Phys. Chem. 2017, 36, 287–331. [CrossRef] 6. Ceccarelli, C.; Caselli, P.; Fontani, F.; Neri, R.; López-Sepulcre, A.; Codella, C.; Feng, S.; Jiménez-Serra, I.; Lefloch, B.; Pineda, J.E.; et al. Seeds Of Life In Space (SOLIS): The Organic Composition Diversity at 300–1000 au Scale in Solar-type Star-forming Regions. Astrophys. J. 2017, 850, 176. [CrossRef] 7. Caselli, P.; Ceccarelli, C. Our astrochemical heritage. Astron. Astrophys. Rev. 2012, 20, 56. [CrossRef] 8. Rotelli, L.; Trigo-Rodríguez, J.M.; Moyano-Cambero, C.E.; Carota, E.; Botta, L.; Di Mauro, E.; Saladino, R. The key role of meteorites in the formation of relevant prebiotic molecules in a formamide/water environment. Sci. Rep. 2016, 6, 38888. [CrossRef] 9. Kwok, S. Complex organics in space from Solar System to distant galaxies. Astron. Astrophys. Rev. 2016, 24, 8. [CrossRef] 10. Kitadai, N.; Maruyama, S. Origins of building blocks of life: A review. Geosci. Front. 2018, 9, 1117–1153. 11. Garrod, R.T.; Herbst, E. Formation of methyl formate and other organic species in the warm-up phase of hot molecular cores. Astron. Astrophys. 2006, 457, 927–936. [CrossRef] 12. Garrod, R.T.; Weaver, S.L.W.; Herbst, E. Complex chemistry in star-forming regions: An expanded gas-grain warm-up chemical model. Astrophys. J. 2008, 682, 283. [CrossRef] 13. Vasyunin, A.I.; Herbst, E. Reactive Desorption and Radiative Association as Possible Drivers of Complex Molecule Formation in the Cold Interstellar Medium. Astrophys. J. 2013, 769, 34. [CrossRef] 14. Balucani, N.; Ceccarelli, C.; Taquet, V. Formation of complex organic molecules in cold objects: The role of gas-phase reactions. Mon. Not. R. Astron. Soc. 2015, 449, L16–L20. [CrossRef] 15. Jørgensen, J.K.; Belloche, A.; Garrod, R.T. Astrochemistry during the formation of stars. Ann. Rev. Astron. Astrophys. 2020, 58, 727–778. [CrossRef] 16. Ceccarelli, C.; Codella, C.; Balucani, N.; Bockelee-Morvan, D.; Herbst, E.; Vastel, C.; Caselli, P.; Favre, C.; Lefloch, B.; Oberg, K.; et al. Organic Chemistry in the First Phases of Solar-Type Protostars. In Proceedings of the Astronomical Society of the Pacific Conference Series; Inutsuka, S., Aikawa, Y., Muto, T., Tomida, K., Tamura, M., Eds.; Cornell University: Ithaca, NY, USA, 2023; Volume 534, p. 379. 17. Tielens, A.G.G.M.; Hagen, W. Model calculations of the molecular composition of interstellar grain mantles. Astron. Astrophys. 1982, 114, 245–260. 18. Molpeceres, G.; Rimola, A.; Ceccarelli, C.; Kästner, J.; Ugliengo, P.; Maté, B. Silicate-mediated interstellar water formation: A theoretical study. Mon. Not. R. Astron. Soc. 2019, 482, 5389–5400. [CrossRef] 19. Ferrero, S.; Pantaleone, S.; Ceccarelli, C.; Ugliengo, P.; Sodupe, M.; Rimola, A. Where does the energy go during the interstellar NH3 formation on water ice? A computational study. Astrophys. J. 2023, 944, 142. [CrossRef] 20. Rimola, A.; Taquet, V.; Ugliengo, P.; Balucani, N.; Ceccarelli, C. Combined quantum chemical and modeling study of CO hydrogenation on water ice. Astron. Astrophys. 2014, 572, A70. [CrossRef] 21. Skouteris, D.; Balucani, N.; Ceccarelli, C.; Vazart, F.; Puzzarini, C.; Barone, V.; Codella, C.; Lefloch, B. The genealogical tree of ethanol: Gas-phase formation of glycolaldehyde, acetic acid, and formic acid. Astrophys. J. 2018, 854, 135. [CrossRef] Int. J. Mol. Sci. 2023, 24, 16824 22 of 25 22. Vazart, F.; Ceccarelli, C.; Balucani, N.; Bianchi, E.; Skouteris, D. Gas-phase formation of acetaldehyde: Review and new theoretical computations. Mon. Not. R. Astron. Soc. 2020, 499, 5547–5561. [CrossRef] 23. Heard, D.E. Rapid acceleration of hydrogen atom abstraction reactions of OH at very low temperatures through weakly bound complexes and tunneling. Acc. Chem. Res. 2018, 51, 2620–2627. [CrossRef] [PubMed] 24. Barone, V.; Latouche, C.; Skouteris, D.; Vazart, F.; Balucani, N.; Ceccarelli, C.; Lefloch, B. Gas-phase formation of the prebiotic molecule formamide: Insights from new quantum computations. Mon. Not. R. Astron. Soc. 2015, 453, L31–L35. [CrossRef] 25. Skouteris, D.; Vazart, F.; Ceccarelli, C.; Balucani, N.; Puzzarini, C.; Barone, V. New quantum chemical computations of formamide deuteration support gas-phase formation of this prebiotic molecule. Mon. Not. R. Astron. Soc. 2017, 468, L1–L5. [CrossRef] 26. Enrique-Romero, J.; Rimola, A.; Ceccarelli, C.; Balucani, N. The (impossible?) formation of acetaldehyde on the grain surfaces: Insights from quantum chemical calculations. Mon. Not. R. Astron. Soc. 2016, 459, L6–L10. [CrossRef] 27. Rimola, A.; Skouteris, D.; Balucani, N.; Ceccarelli, C.; Enrique-Romero, J.; Taquet, V.; Ugliengo, P. Can Formamide Be Formed on In- 697 terstellar Ice? An Atomistic Perspective. ACS Earth Space Chem. 2018, 2, 720–734. [CrossRef] 28. Enrique-Romero, J.; Ceccarelli, C.; Rimola, A.; Skouteris, D.; Balucani, N.; Ugliengo, P. Theoretical computations on the efficiency of acetaldehyde formation on interstellar icy grains. Astron. Astrophys. 2021, 655, A9. [CrossRef] 29. Enrique-Romero, J.; Rimola, A.; Ceccarelli, C.; Ugliengo, P.; Balucani, N.; Skouteris, D. Quantum Mechanical Simulations of the Radical–Radical Chemistry on Icy Surfaces. Astrophys. J. Suppl. Ser. 2022, 259, 39. [CrossRef] 30. Molpeceres, G.; Jiménez-Serra, I.; Oba, Y.; Nguyen, T.; Watanabe, N.; de la Concepción, J.G.; Maté, B.; Oliveira, R.; Kästner, J. Hydrogen abstraction reactions in formic and thioformic acid isomers by hydrogen and deuterium atoms. Astron. Astrophys. 2022, 663, A41. [CrossRef] 31. Jin, M.; Garrod, R.T. Formation of complex organic molecules in cold interstellar environments through nondiffusive grain-surface and ice-mantle chemistry. Astrophys. J. Suppl. Ser. 2020, 249, 26. [CrossRef] 32. Garrod, R.T.; Jin, M.; Matis, K.A.; Jones, D.; Willis, E.R.; Herbst, E. Formation of complex organic molecules in hot molecular cores through nondiffusive grain-surface and ice-mantle chemistry. Astrophys. J. Suppl. Ser. 2022, 259, 1. [CrossRef] 33. Molpeceres, G.; Kastner, J.; Fedoseev, G.; Qasim, D.; Schomig, R.; Linnartz, H.; Lamberts, T. Carbon atom reactivity with amorphous solid water: H2 O-catalyzed formation of H2 CO. J. Phys. Chem. Lett. 2021, 12, 10854–10860. [CrossRef] [PubMed] 34. Perrero, J.; Enrique-Romero, J.; Martínez-Bachs, B.; Ceccarelli, C.; Balucani, N.; Ugliengo, P.; Rimola, A. Non-energetic formation of ethanol via CCH reaction with interstellar H2 O ices. A computational chemistry study. ACS Earth Space Chem. 2022, 6, 496–511. [CrossRef] [PubMed] 35. Perrero, J.; Ugliengo, P.; Ceccarelli, C.; Rimola, A. Quantum mechanical modelling of the grain-surface formation of acetaldehyde on H2 O: CO dirty ice surfaces. Mon. Not. R. Astron. Soc. 2023, 525, 2654–2667. [CrossRef] 36. Rubin, R.; Swenson, G., Jr.; Benson, R.; Tigelaar, H.; Flygare, W. Microwave detection of interstellar formamide. Astrophys. J. 1971, 169, L39. [CrossRef] 37. Kahane, C.; Ceccarelli, C.; Faure, A.; Caux, E. Detection of formamide, the simplest but crucial amide, in a solar-type protostar. Astrophys. J. Lett. 2013, 763, L38. [CrossRef] 38. Mendoza, E.; Lefloch, B.; López-Sepulcre, A.; Ceccarelli, C.; Codella, C.; Boechat-Roberty, H.; Bachiller, R. Molecules with a peptide link in protostellar shocks: A comprehensive study of L1157. Mon. Not. R. Astron. Soc. 2014, 445, 151–161. [CrossRef] 39. Jaber, A.A.; Ceccarelli, C.; Kahane, C.; Caux, E. The census of complex organic molecules in the solar-type protostar IRAS16293- 2422. Astrophys. J. 2014, 791, 29. [CrossRef] 40. López-Sepulcre, A.; Jaber, A.A.; Mendoza, E.; Lefloch, B.; Ceccarelli, C.; Vastel, C.; Bachiller, R.; Cernicharo, J.; Codella, C.; Kahane, C.; et al. Shedding light on the formation of the pre-biotic molecule formamide with ASAI. Mon. Not. R. Astron. Soc. 2015, 449, 2438–2458. [CrossRef] 41. Marcelino, N.; Gerin, M.; Cernicharo, J.; Fuente, A.; Wootten, H.; Chapillon, E.; Pety, J.; Lis, D.C.; Roueff, E.; Commerçon, B.; et al. ALMA observations of the young protostellar system Barnard 1b: Signatures of an incipient hot corino in B1b-S. Astron. Astrophys. 2018, 620, A80. [CrossRef] 42. Bianchi, E.; Codella, C.; Ceccarelli, C.; Vazart, F.; Bachiller, R.; Balucani, N.; Bouvier, M.; De Simone, M.; Enrique-Romero, J.; Kahane, C.; et al. The census of interstellar complex organic molecules in the Class I hot corino of SVS13-A. Mon. Not. R. Astron. Soc. 2019, 483, 1850–1861. [CrossRef] 43. López-Sepulcre, A.; Balucani, N.; Ceccarelli, C.; Codella, C.; Dulieu, F.; Theulé, P. Interstellar formamide (NH2 CHO), a key prebiotic precursor. ACS Earth Space Chem. 2019, 3, 2122–2137. [CrossRef] 44. Song, L.; Kästner, J. Formation of the prebiotic molecule NH2 CHO on astronomical amorphous solid water surfaces: Accurate tunneling rate calculations. Phys. Chem. Chem. Phys. 2016, 18, 29278–29285. [CrossRef] [PubMed] 45. Enrique-Romero, J.; Rimola, A.; Ceccarelli, C.; Ugliengo, P.; Balucani, N.; Skouteris, D. Reactivity of HCO with CH3 and NH2 on water ice surfaces. a comprehensive accurate quantum chemistry study. ACS Earth Space Chem. 2019, 3, 2158–2170. [CrossRef] 46. Vazart, F.; Calderini, D.; Puzzarini, C.; Skouteris, D.; Barone, V. State-of-the-art thermochemical and kinetic computations for astrochemical complex organic molecules: Formamide formation in cold interstellar clouds as a case study. J. Chem. Theory Comput. 2016, 12, 5385–5397. [CrossRef] [PubMed] 47. Gottlieb, C.A. Detection of acetaldehyde in Sagittarius. In Molecules in the Galactic Environment; John Wiley and Sons: Hoboken, NJ, USA, 1973; p. 181. 48. Pizzarello, S.; Weber, A.L. Prebiotic amino acids as asymmetric catalysts. Science 2004, 303, 1151. [CrossRef] [PubMed] Int. J. Mol. Sci. 2023, 24, 16824 23 of 25 49. Córdova, A.; Ibrahem, I.; Casas, J.; Sundén, H.; Engqvist, M.; Reyes, E. Amino acid catalyzed neogenesis of carbohydrates: A plausible ancient transformation. Chem. Eur. J. 2005, 11, 4772–4784. [CrossRef] 50. Charnley, S. Acetaldehyde in star-forming regions. Adv. Space Res. 2004, 33, 23–30. [CrossRef] 51. Vasyunin, A.I.; Caselli, P.; Dulieu, F.; Jiménez-Serra, I. Formation of complex molecules in prestellar cores: A multilayer approach. Astrophys. J. 2017, 842, 33. [CrossRef] 52. Lamberts, T.; Markmeyer, M.N.; Kolb, F.J.; Kastner, J. Formation of acetaldehyde on CO-rich ices. ACS Earth Space Chem. 2019, 3, 958–963. [CrossRef] 53. Enrique-Romero, J.; Álvarez-Barcia, S.; Kolb, F.; Rimola, A.; Ceccarelli, C.; Balucani, N.; Meisner, J.; Ugliengo, P.; Lamberts, T.; Kästner, J. Revisiting the reactivity between HCO and CH3 on interstellar grain surfaces. Mon. Not. R. Astron. Soc. 2020, 493, 2523–2527. [CrossRef] 54. Ferrero, S.; Ceccarelli, C.; Ugliengo, P.; Sodupe, M.; Rimola, A. Formation of complex organic molecules on interstellar CO ices? Insights from computational chemistry simulations. Astrophys. J. 2023, 951, 150. [CrossRef] 55. Xu, S.; Lin, M.C. Theoretical study on the kinetics for OH reactions with CH3 OH and C2 H5 OH. Proc. Combust. Inst. 2007, 31, 159–166. [CrossRef] 56. Galano, A.; Alvarez-Idaboy, J.R.; Bravo-Pérez, G.; Ruiz-Santoyo, M.E. Gas phase reactions of C1–C4 alcohols with the OH radical: A quantum mechanical approach. Phys. Chem. Chem. Phys. 2002, 4, 4648–4662. [CrossRef] 57. Zheng, J.; Truhlar, D.G. Multi-path variational transition state theory for chemical reaction rates of complex polyatomic species: Ethanol + OH reactions. Faraday Discuss. 2012, 157, 59–88. [CrossRef] [PubMed] 58. Marinov, N.M. A detailed chemical kinetic model for high temperature ethanol oxidation. Int. J. Chem. Kinet. 1999, 31, 183–220. [CrossRef] 59. Carr, S.A.; Blitz, M.A.; Seakins, P.W. Site-specific rate coefficients for reaction of OH with ethanol from 298 to 900 K. J. Phys. Chem. A 2011, 115, 3335–3345. [CrossRef] [PubMed] 60. Brown, R.; Crofts, J.; Gardner, F.; Godfrey, P.; Robinson, B.; Whiteoak, J. Discovery of interstellar methyl formate. Astrophys. J. 1975, 197, L29–L31. [CrossRef] 61. Hollis, J.; Vogel, S.; Snyder, L.; Jewell, P.; Lovas, F. The spatial scale of glycolaldehyde in the galactic center. Astrophys. J. 2001, 554, L81. [CrossRef] 62. Horn, A.; Møllendal, H.; Sekiguchi, O.; Uggerud, E.; Roberts, H.; Herbst, E.; Viggiano, A.; Fridgen, T.D. The gas-phase formation of methyl formate in hot molecular cores. Astrophys. J. 2004, 611, 605. [CrossRef] 63. Laas, J.C.; Garrod, R.T.; Herbst, E.; Weaver, S.L.W. Contributions from grain surface and gas phase chemistry to the formation of methyl formate and its structural isomers. Astrophys. J. 2011, 728, 71. [CrossRef] 64. Ruaud, M.; Loison, J.; Hickson, K.C.; Gratier, P.M.; Hersant, F.W.V. Modelling complex organic molecules in dense regions: Eley-Rideal and complex induced reaction. Mon. Not. R. Astron. Soc. 2015, 447, 4004–4017. [CrossRef] 65. Hoyermann, K.; Nacke, F. Elementary reactions of the methoxymethyl radical in the gas phase: CH3 OCH3 +F, CH2 OCH3 +CH2 OCH3 , CH2 OCH3 +O2 and CH2 OCH3 +O. Symp. Int. Combust. Proc. 1996, 26, 505–512. [CrossRef] 66. Song, X.; Hou, H.; Wang, B. Mechanistic and kinetic study of the O + CH3 OCH2 reaction and the unimolecular decomposition of CH3 OCH2 O. Phys. Chem. Chem. Phys. 2005, 7, 3980–3988. [CrossRef] [PubMed] 67. Zuckerman, B.; Ball, J.A.; Gottlieb, C.A. Microwave detection of interstellar formic acid. Astrophys. J. 1971, 163, L41. [CrossRef] 68. Pilling, S.; Baptista, L.; Boechat-Roberty, H.M.; Andrade, D.P. Formation routes of interstellar glycine involving carboxylic acids: Possible favoritism between gas and solid phase. Astrobiology 2011, 11, 883–893. [CrossRef] [PubMed] 69. Schmidt, F.; Swiderek, P.; Bredehoft, J.H. Formation of formic acid, formaldehyde, and carbon dioxide by electron-induced chemistry in ices of water and carbon monoxide. ACS Earth Space Chem. 2019, 3, 1974–1986. [CrossRef] 70. Woon, D.E. Quantum Chemical Cluster Studies of Cation–Ice Reactions for Astrochemical Applications: Seeking Experimental Confirmation. Acc. Chem. Res. 2021, 54, 490–497. [CrossRef] 71. McElroy, D.; Walsh, C.; Markwick, A.; Cordiner, M.; Smith, K.; Millar, T. The UMIST database for astrochemistry 2012. Astron. Astrophys. 2013, 550, A36. [CrossRef] 72. Wakelam, V.; Loison, J.C.; Herbst, E.; Pavone, B.; Bergeat, A.; Béroff, K.; Chabot, M.; Faure, A.; Galli, D.; Geppert, W.D.; et al. The 2014 KIDA network for interstellar chemistry. Astrophys. J. Suppl. Ser. 2015, 217, 20. [CrossRef] 73. Yetter, R.A.; Rabitz, H.; Dryer, F.L.; Maki, R.G.; Klemm, R.B. Evaluation of the rate constant for the reaction OH+H2 CO: Application of modeling and sensitivity analysis techniques for determination of the product branching ratio. J. Chem. Phys. 1989, 91, 4088–4097. [CrossRef] 74. D’Anna, B.; Bakken, V.; Beukes, J.A.; Nielsen, C.J.; Brudnik, K.; Jodkowski, J.T. Experimental and theoretical studies of gas phase NO3 and OH radical reactions with formaldehyde, acetaldehyde and their isotopomers. Phys. Chem. Chem. Phys. 2003, 5, 1790–1805. [CrossRef] 75. Xu, S.; Zhu, R.; Lin, M.C. Ab initio study of the OH+ CH2 O reaction: The effect of the OH···OCH2 complex on the H-abstraction kinetics. Int. J. Chem. Kinet. 2006, 38, 322–326. [CrossRef] 76. Garrod, R.T. A three-phase chemical model of hot cores: The formation of glycine. Astrophys. J. 2013, 765, 60. [CrossRef] 77. Collings, M.P.; Frankland, V.L.; Lasne, J.; Marchione, D.; Rosu-Finsen, A.; McCoustra, M.R.S. Probing model interstellar grain surfaces with small molecules. Mon. Not. R. Astron. Soc. 2015, 449, 1826–1833. [CrossRef] Int. J. Mol. Sci. 2023, 24, 16824 24 of 25 78. Yu, A.; Tatiana, F.; Drebushchak, N.; Tantardini, C. Seeking the best model for non-covalent interactions within the crystal structure of meloxicam. Comput. Theor. Chem. 2019, 1157, 47–53. [CrossRef] 79. Varadwaj, P.R.; Varadwaj, A.; Marques, H.M.; Yamashita, K. Significance of hydrogen bonding and other noncovalent interactions in determining octahedral tilting in the CH3 NH3 PbI3 hybrid organic-inorganic halide perovskite solar cell semiconductor. Sci. Rep. 2019, 9, 50. [CrossRef] 80. Chaudhari, S.; Nikam, S.; Khatri, N.; Wakde, S. Co-crystals: A review. J. Drug Deliv. Ther. 2018, 8, 350–358. [CrossRef] 81. Fucke, K.; Myz, S.A.; Shakhtshneider, T.P.; Boldyreva, E.V.; Griesser, U.J. How good are the crystallisation methods for co-crystals? A comparative study of piroxicam. New J. Chem. 2012, 36, 1969–1977. [CrossRef] 82. Minissale, M.; Aikawa, Y.; Bergin, E.; Bertin, M.; Brown, W.A.; Cazaux, S.; Charnley, S.B.; Coutens, A.; Cuppen, H.M.; Guzman, V.; et al. Thermal desorption of interstellar ices: A review on the controlling parameters and their implications from snowlines to chemical complexity. ACS Earth Space Chem. 2022, 6, 597–630. [CrossRef] 83. Congiu, E.; Minissale, M.; Baouche, S.; Chaabouni, H.; Moudens, A.; Cazaux, S.; Manicò, G.; Pirronello, V.; Dulieu, F. Efficient diffusive mechanisms of O atoms at very low temperatures on surfaces of astrophysical interest. Faraday Discuss. 2014, 168, 151. [CrossRef] [PubMed] 84. Lee, M.W.; Meuwly, M. Diffusion of atomic oxygen relevant to water formation in amorphous interstellar ices. Faraday Discuss. 2014, 168, 205. [CrossRef] [PubMed] 85. Perrero, J.; Enrique-Romero, J.; Ferrero, S.; Ceccarelli, C.; Podio, L.; Codella, C.; Rimola, A.; Ugliengo, P. Binding Energies of Interstellar Relevant S-bearing Species on Water Ice Mantles: A Quantum Mechanical Investigation. Astrophys. J. 2022, 938, 158. [CrossRef] 86. Pantaleone, S.; Enrique-Romero, J.; Ceccarelli, C.; Ugliengo, P.; Balucani, N.; Rimola, A. Chemical desorption versus energy dissipation: Insights from ab initio molecular dynamics of HCO· formation. Astrophys. J. 2020, 897, 56. [CrossRef] 87. Pantaleone, S.; Enrique-Romero, J.; Ceccarelli, C.; Ferrero, S.; Balucani, N.; Rimola, A.; Ugliengo, P. H2 formation on interstellar grains and the fate of reaction energy. Astrophys. J. 2021, 917, 49. [CrossRef] 88. Balakrishnan, N. Quantum mechanical investigation of the O + H2 → OH + H reaction. J. Chem. Phys. 2003, 119, 195–199. [CrossRef] 89. Meisner, J.; Kamp, I.; Thi, W.F.; Kästner, J. The role of atom tunneling in gas-phase reactions in planet-forming disks. Astron. Astrophys. 2019, 627, A45. [CrossRef] 90. Fermann, J.T.; Auerbach, S. Modeling proton mobility in acidic zeolite clusters: II. Room temperature tunneling effects from semiclassical rate theory. J. Chem. Phys. 2000, 112, 6787–6794. [CrossRef] 91. Bao, J.L.; Truhlar, D.G. Variational transition state theory: Theoretical framework and recent developments. Chem. Soc. Rev. 2017, 46, 7548–7596. [CrossRef] 92. Kästner, J. Theory and simulation of atom tunneling in chemical reactions. WIREs Comput. Mol. Sci. 2014, 4, 158–168. [CrossRef] 93. Pisani, C.; Casassa, S.; Ugliengo, P. Proton-ordered ice structures at zero pressure. A quantum-mechanical investigation. Chem. Phys. Lett. 1996, 253, 201–208. [CrossRef] 94. Casassa, S.; Ugliengo, P.; Pisani, C. Proton-ordered models of ordinary ice for quantum-mechanical studies. J. Chem. Phys. 1997, 106, 8030–8040. [CrossRef] 95. Zamirri, L.; Casassa, S.; Rimola, A.; Segado-Centellas, M.; Ceccarelli, C.; Ugliengo, P. IR spectral fingerprint of carbon monoxide in interstellar water–ice models. Mon. Not. R. Astron. Soc. 2018, 480, 1427–1444. [CrossRef] 96. Watanabe, N.; Kouchi, A. Ice surface reactions: A key to chemical evolution in space. Prog. Surf. Sci. 2008, 83, 439–489. [CrossRef] 97. Ferrero, S.; Zamirri, L.; Ceccarelli, C.; Witzel, A.; Rimola, A.; Ugliengo, P. Binding Energies of Interstellar Molecules on Crystalline and Amorphous Models of Water Ice by Ab Initio Calculations. Astrophys. J. 2020, 904, 11. [CrossRef] 98. Shimonishi, T.; Nakatani, N.; Furuya, K.; Hama, T. Adsorption Energies of Carbon, Nitrogen, and Oxygen Atoms on the Low-temperature Amorphous Water Ice: A Systematic Estimation from Quantum Chemistry Calculations. Astrophys. J. 2018, 732, 27. [CrossRef] 99. Dovesi, R.; Erba, A.; Orlando, R.; Zicovich-Wilson, C.M.; Civalleri, B.; Maschio, L.; Rérat, M.; Casassa, S.; Baima, J.; Salustro, S.; et al. Quantum-mechanical condensed matter simulations with CRYSTAL. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2018, 8, e1360. [CrossRef] 100. Sure, R.; Grimme, S. Corrected small basis set Hartree-Fock method for large systems. J. Comput. Chem. 2013, 34, 1672. [CrossRef] [PubMed] 101. Kruse, H.; Grimme, S. A geometrical correction for the inter-and intra-molecular basis set superposition error in Hartree-Fock and density functional theory calculations for large systems. J. Chem. Phys. 2012, 136, 154101. [CrossRef] 102. Tatewaki, H.; Huzinaga, S. A systematic preparation of new contracted Gaussian-type orbital sets. III. Second-row atoms from Li through Ne. J. Comput. Chem. 1980, 1, 205. [CrossRef] 103. Liu, B.; McLean, A. Accurate calculation of the attractive interaction of two ground state helium atoms. J. Chem. Phys. 1973, 59, 4557. [CrossRef] 104. Jansen, H.; Ros, P. Non-empirical molecular orbital calculations of the protonation of carbon monoxide. Chem. Phys. Lett. 1969, 3, 140. [CrossRef] 105. Becke, A.D. A new mixing of Hartree–Fock and local density-functional theories. J. Chem. Phys. 1993, 98, 1372–1377. [CrossRef] 106. Lee, C.; Yang, W.; Parr, R.G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 1988, 37, 785–789. [CrossRef] Int. J. Mol. Sci. 2023, 24, 16824 25 of 25 107. Zhao, Y.; Truhlar, D.G. The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncova- lent interactions, excited states, and transition elements: Two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor. Chem. Acc. 2008, 120, 215–241. 108. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. [CrossRef] 109. Grimme, S.; Ehrlich, S.; Goerigk, L. Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 2011, 32, 1456–1465. [CrossRef] 110. Frisch, M.J.; Trucks, G.W.; Schlegel, H.B.; Scuseria, G.E.; Robb, M.A.; Cheeseman, J.R.; Scalmani, G.; Barone, V.; Petersson, G.A.; Nakatsuji, H.; et al. Gaussian 16, Revision C.01; Gaussian Inc.: Wallingford, CT, USA, 2016. 111. Dunning, T.H., Jr. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. J. Chem. Phys. 1989, 90, 1007–1023. [CrossRef] 112. Pascale, F.; Zicovich-Wilson, C.M.; López Gejo, F.; Civalleri, B.; Orlando, R.; Dovesi, R. The calculation of the vibrational frequencies of crystalline compounds and its implementation in the CRYSTAL code. J. Comput. Chem. 2004, 25, 888–897. [CrossRef] 113. Zicovich-Wilson, C.M.; Pascale, F.; Roetti, C.; Saunders, V.R.; Orlando, R.; Dovesi, R. Calculation of the vibration frequencies of quartz: The effect of Hamiltonian and basis set. J. Comput. Chem. 2004, 25, 1873–1881. [CrossRef] Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.