Theoretical Determination of Binding Energies of Small Molecules on Interstellar Ice Surfaces

The adsorption of a series of atoms and small molecules and radicals (H, C, N, O, NH, OH, H2O, CH3, and NH3) on hexagonal crystalline and amorphous ice clusters were obtained via classical molecular dynamics and electronic structure methods. The geometry and binding energies were calculated using a QMHigh:QMLow hybrid method on model clusters. Several combination of basis sets, density functionals and semi-empirical methods were compared and tested against previous works. More accurate binding energies were also refined via single point Coupled Cluster calculations. Most species, except carbon atom, physisorb on the surface, leading to rather small binding energies. The carbon atom forms a COH2 molecule and in some cases leads to the formation of a COH-H3O+ complex. Amorphous ices are characterized by slightly stronger binding energies than the crystalline phase. A major result of this work is to also access the dispersion of the binding energies since a variety of adsorption sites is explored. The interaction energies thus obtained may serve to feed or refine astrochemical models. The present methodology could be easily extended to other types of surfaces and larger adsorbates.


INTRODUCTION
Astrochemistry has been an area of intense research in particular with the aim to elucidate the origin of the chemical diversity observed in space. Depending on the evolution stage of the interstellar dust and the region, energetic processes may dominate the chemistry. However, in low UV-photon field regions, as in dense molecular clouds, energetic processes are suppressed and surface reactions of cold molecules become important. It has been demonstrated both experimentally and theoretically that dust grains catalyze the formation of some important molecules, such as H 2 , H 2 O, and CO 2 , organic molecules (H 2 CO, CH 3 OH . . . ) (Ioppolo et al., 2011;Pirim and Krim, 2011;Peters et al., 2013;Rimola et al., 2014;Minissale et al., 2015), or complex/prebiotic molecules (Koch et al., 2008;Rimola and Ugliengo, 2009). Thermally activated surface reactions rarely occur over the activation energy barriers on ice mantles because of their low temperatures (∼ 10 K). In contrast, quantumtunneling reactions are facilitated in this temperature range. The importance of tunneling effects has been discussed in various studies (Goumans and Kästner, 2010;Hama and Watanabe, 2013;Hickson et al., 2016). Whatever the conditions are, for a reaction to occur on the grain, the molecular species needs to reside or diffuse at the grain surface. Most neutral species can physisorb onto grains through van der Waals interactions at the low temperature conditions of dense molecular clouds. Chemisorption may also occur for radicals or reactive molecules on specific or defective sites. In contrast to chemisorption, physisorption is reversible, the desorption rate being derived from the binding energy and the frequency characterizing the motion of the adsorbed species in the potential energy well. Diffusion is also to be considered when reactions or recombinations occur through Langmuir-Hinshelwood mechanisms. Most astrochemical models Taquet et al., 2013;Ruaud et al., 2016;Cuppen et al., 2017;Wakelam et al., 2017) mimicking the chemical evolution in star forming region have now integrated these gas-grain processes. Among the molecular data used as inputs to these models, the adsorption energy is a key parameter, the diffusion energy being estimated as a fraction (between 0.3 and 1.0) of the binding energy (Ruaud et al., 2016). There is thus a crucial need to improve the binding energies (BE's) values in the models to better reproduce the astrophysical abundances.
Binding energies depend indeed on the nature of the volatile species, the surface and its coverage. BE's can be determined experimentally for stable molecules using programmed desorption (TPD) methods (Burke and Brown, 2010;Hama and Watanabe, 2013). Several studies have reviewed known binding energies on various surfaces and more especially water ices (Fraser et al., 2001;Collings et al., 2004;Noble et al., 2012;Hama and Watanabe, 2013;Penteado et al., 2017;Wakelam et al., 2017). However, TPD is not appropriate for reactive atoms or radicals, since they may spontaneously recombine to form stable molecules at the surface before thermal desorption. Consequently, binding energies of atoms or radicals are less documented.
For species with undefined BE's, attempts have been made by speculating an additive law: the BE of a given molecule being the sum of the BE's of the atoms or groups of atoms. This approximation of the BE as implemented in Garrod et al.'s model (Cuppen et al., 2017) gives a rough estimation that can indeed be improved in different manners. Wakelam et al. (Wakelam et al., 2017) have proposed a proportional law that relates the interaction energy between the given species and one water molecule (energy of the dimer) and its binding energy on amorphous solid water. However, the integration of the binding energies computed with this method in the gas-grain code Nautilus (Ruaud et al., 2016) did not improve the prediction of the gas phase abundances. A more precise determination of the BEs and inclusion of nonthermal processes have been invoked as possible improvements. Because of their accuracy, electronic structure calculations appear to be the method of choice to improve the BE determination. However, for computational reasons, this approach often reduces the explicit description of the surface to a small number of atoms or molecules (Koch et al., 2007;Duvernay et al., 2014;Rimola et al., 2014). Thanks to the development of numerical resources, periodic ab initio simulations have been also performed recently to account for larger crystalline ice surfaces though still in a static way neglecting temperature effects and the configurations sampling being limited to the identification of the global minimum (Ferrero et al., 2020a). It should be noted that the size of the crystal and the choice of the functional are critical to surely assess an accurate value of the binding energy (Ferrero et al., 2020b). By contrast, semi-empirical based classical molecular dynamics (MD) simulations account at low computational cost for long range effects through the modeling of quasi-infinite surfaces and offer a full exploration of the possible binding sites providing BE distributions and eventually diffusion barriers (Michoulier et al., 2018). However, classical MD rely strongly on the accuracy, or even availability, of force-field parameters. Such parameters are in particular not available for radicals. Hybrid QM/MM (Quantum Mechanics/Molecular Mechanics) approaches, such as the ONIOM method (Chung et al., 2015), recently used by Sameera et al. (Sameera et al., 2017;Sameera and Maseras, 2018) on crystalline water ices, have the advantage of describing accurately the adsorbate-surface interaction including long range effects, dividing the system into two subunits: the adsorbate and few molecules of the surface treated at a high QM level and the rest of the surface being treated with a MM or lower level QM' approach. The latter QM/QM' scheme is the methodology chosen in the present study on ice surfaces generated from classical molecular dynamics simulations, including amorphous ice samples, to quantify the effect of structural disorder.
The binding energy depends on the volatile but also on the grain surface composition. Numerous studies have been conducted to determine the composition of the grains and more especially the grain icy mantle (Whittet et al., 1996;van Dishoeck, 1998;Ehrenfreund and Charnley, 2000;Ehrenfreund et al., 2003;Pontoppidan et al., 2008). They revealed that interstellar ices are mainly composed of water, but may also contain other species like CO, CO 2 , NH 3 , and methanol, in less significant quantities. The precise structure of interstellar ice remains an open question, though it is generally admitted that it might be amorphous (Dulieu et al., 2013;Hama and Watanabe, 2013). Amorphous Solid Water (ASW) is known to have physical properties that are distinct from other crystalline states, e.g., a lower thermal conductivity, as well as a larger surface area and higher porosity (Dohnálek et al., 2003). In the spectrum of ASW, a small feature is present on the high wavenumber wing of the OH stretching modes band, due to OH dangling bonds, that are more abundant in the microporous material than in crystalline compact ice (Noble et al., 2014). These unsaturated sites are preferential adsorption sites for incoming molecules (Buch and Devlin, 1991;Michoulier et al., 2020). Proton disorder in amorphous ice was shown to trigger a diversity of adsorption sites and consequently of binding energies for small adsorbates compared to crystalline ice (Kimmel et al., 2001;Kimmel et al., 2001). Even though laboratory ASW or ice grown in MD simulations do not reproduce interstellar icy grains for various factors (growth conditions, porosity, size and time scales . . . ) (Hama and Watanabe, 2013), characterizing the BEs of various molecules on these model surfaces provide nonetheless valuable insights on the trapping of the gas-phase species, on their availability for surface reactions and on the influence of the surface heterogeneity on the BE values.
The purpose of the present study is to establish a theoretical framework combining a force field based MD sampling and density functional corrected with high level ab initio calculations on some configurations issued from this sampling to accurately describe binding energies (physisorption or chemisorption) of a set of atoms, radicals or molecules on surfaces of interstellar interest. This consists in probing a rich variability of adsorption sites on crystalline or amorphous ices, each of which characterized by a specific binding energy distribution.
A first section presents the details of the methodology and of the performed calculations. Results are then presented. The crystalline ice investigation serves firstly as a benchmark of the model and QM method. Once set, the most accurate combination of methods is then applied on amorphous ices. The obtained values are discussed and compared to published data when available.

METHODS
As a start, it is important to agree on the definition of the BE of a species on a substrate. Throughout the paper, BE will be defined as the positive quantity according to the following equation: Where E ads is the total calculated energy of the species adsorbed on the surface while E substrate and E species are obtained separately. Within the goal to describe adsorption on both periodic and amorphous ices, the ONIOM hybrid method has been chosen (Chung et al., 2015) as implemented in the Gaussian 16 code (Frisch et al., 2016). To avoid problems inherent to QM/MM methods, especially in the Mechanical Embedding scheme previously employed (Sameera et al., 2017;Sameera and Maseras, 2018;Sameera et al., 2021), we rather chose a QMHigh:QMLow combination. The QMHigh part was described with Density Functional Theory (DFT). Following previous works from Sameera et al. (Sameera et al., 2017; FIGURE 1 | Top view of a snapshot of amorphous ice issued from a MD trajectory at 77 K. The layers for the ONIOM calculation are represented accordingly: adsorbed H 2 O (beige), QMHigh (red/white), QMLow (blue), and the remaining water molecules (grey). Sameera and Maseras, 2018), two gaussian atomic basis sets were used, namely Pople's 6-31+G** (Ditchfield et al., 1971;Hehre et al., 1972;Hariharan and Pople, 1973;Clark et al., 1983) and the def2-TZVP (Weigend and Ahlrichs, 2005). The same functionals as in Sameera et al. were also employed: the M06-2X (Zhao and Truhlar, 2008) and the range separated ωB97X-D (Chai and Head-Gordon, 2008), with the Gaussian 16 default grid size (i.e., ultrafine). For the QMLow part, we used PM6 (Stewart, 2007) and PM7R8 (Throssell and Frisch, 2018;Throssell, 2017) parameters. According to Throssel (Throssell, 2017), the PM7R8 parametrization improves the description of hydrogen bonded systems. For open shell systems, i.e., atoms and radicals, only the high-spin ground state was considered and the unrestricted formalism was employed.
Geometry optimization were carried out in Cartesian coordinates on the QMHigh part, the QMLow zone remaining frozen. Default convergence thresholds on gradients and displacements were used. It should be noticed that no BSSE effect was considered in the present work. For crystalline (Ih) ices, the geometries were obtained from the supplementary material of Sameera et al. (Sameera et al., 2017). They are derived from Karseemeijer et al. (Karssemeijer et al., 2012) and were obtained at the ONIOM(DFT:AMOEBA) level. Thus, two ice models with 4 sites differing by the surface dangling bonds ordering were used. In model A, 162 water molecules are present, with 24 in the QMHigh region. Model B comprises 165 water molecules, the QMHigh zone containing 20 molecules (see Figure 2 in Sameera et al. for a graphical representation). For OH, HCO, and CH 3 , the Sameera et al. geometries were used as starting guesses. For the 6 other adsorbates, the initial geometry was generated manually. The harmonic frequencies were obtained at the ωB97X-D/6-31+G**:PM6 level only. The resulting total zero-point vibrational energy correction (ZPE) was then scaled using a factor of 0.9523 taken from the NIST database (https://cccbdb. nist.gov/vsfx.asp). This ZPE correction was then added to the pure electronic BE's obtained will all the different levels of theory employed: For amorphous cases, a different strategy was adopted. An amorphous ice slab was generated by classical molecular dynamics simulations following the procedure detailed in Michoulier et al. (Michoulier et al., 2018), using the TIP4P/ 2005 water model (Abascal and Vega, 2005;Vega and Abascal, 2011). Classical molecular dynamics trajectories have been used to model the adsorption of a single water molecule, described as well with the TIP4P/2005 model, on the generated ice slab at T 77 K. The goal of these independent trajectories is to sample a variety of adsorption sites that will serve as guess sites for the other adsorbates. The MD simulation details can be found in (Michoulier et al., 2018). The amorphous ice slab contains here 1344 molecules in a simulation box of dimensions: 25.89 × 29.95 × 90.28 A. The initial position of the adsorbed water molecule is randomly chosen for each trajectory. After a 1 ns trajectory with a 1 fs timestep at T 77 K, the final configuration of each trajectory is collected, the conservation of total energy and the convergence of the usual parameters being carefully checked for each trajectory.
The obtained adsorption geometries were cut into 3 different layers. The first one defining the QMHigh part, consists in the absorbed molecule surrounded by a sphere containing about 20 other water molecules. Then, about 120-150 molecules in a slightly larger sphere were selected for the QMLow part and finally the remaining ∼1200 ones were removed. The radii of the spheres were adjusted to obtain approximately the same number of molecules as in the crystalline case. The typical radius of these spheres is about 7-8 Å for the High-level part and about 12-15 Å for the Low level one. The reduction of the number of water molecules from the slab to a reduced cluster is illustrated in Figure 1. The same ONIOM procedure used for the crystalline phase was applied to these clusters. For reasons explained below, these calculations were performed at the ωB97X-D/6-31+G**: PM6 level only. The other adsorption cases were treated by replacing the water molecule by the corresponding adsorbate  and reoptimizing the geometries. This method proved to give reasonable guess structures for the ONIOM geometry optimization, despite the different nature of the interaction between the various species and the surface. As in the crystalline case, the QMHigh part without the adsorbate (bare ice surface) was optimized separately as well.
The procedure described in the previous paragraphs gives geometry and binding energies at the DFT level. In order to obtain more accurate binding energies, the ωB97X-D/6-31+G**: PM6 QMHigh partial geometries were extracted and used as inputs for a CBS/DLPNO-CCSD(T) (Guo et al., 2018) single point calculation, using the ORCA code (Neese, 2018;Neese et al., 2020). Resolution of the identity (RIJK) was used (Weigend et al., 2009) with the auxiliary basis set from Weigend (Weigend, 2008). For DLPNO-CCSD(T), the normal criteria were used as discussed by Liakos et al. (Liakos et al., 2015). The auxiliary basis set for  correlation is taken from Hellweg et al. (Hellweg et al., 2007). For Complete Basis Set (CBS) extrapolation, def2-TZVPP and def2-QVZPP basis sets were employed (Weigend and Ahlrichs, 2005). The Hartree-Fock and the correlation energies were extrapolated separately as shown in (Neese and Valeev, 2011) with α 7.880 for Hartree-Fock and β 2.970 for correlation. The three (def2-TZVPP, def2-QVZPP and CBS) resulting DLPNO-CCSD(T) energies then replaced the DFT part of the total ONIOM extrapolated energy. The most sophisticated BE can thus be considered as CBS/DLPNO-CCSD(T):PM6//ωB97X-D/6-31+G**:PM6, corrected with the ωB97X-D/6-31+G**: PM6 ZPE. CBS results should also be free of BSSE (Basis Set Superposition Error), but it should also be noted that BSSE effects were calculated to be small by Sameera et al. (Sameera et al., 2017).
In the following, binding energies are expressed in eV. They can be converted in Kelvins using: 1.000 eV 11604.59 K 96.4853 kJ/mol.

Crystalline ice
The study of the crystalline structures was performed to benchmark the various DFT/ONIOM models employed in the present work. Furthermore, the crystal results may be compared with the amorphous ones as well. We have adopted the nomenclature of Sameera et al. (Sameera et al., 2017) who distinguish 8 different adsorption sites (A1, A2, A3, A4, B1, B2, B3, and B4) depending on the different dangling bonds arrangement at the surface. A detailed comparison with Sameera et al. for the eight sites is given in the supplementary material. In the perspective of studying amorphous ices, we will only report here four characteristic BE's for each species: the minimum, the maximum, the average values and the standard deviation σ, respectively. In the following, we present in the figures the geometries corresponding to the most stable configuration (largest BE)   obtained at ONIOM(ωB97X-D/6-31+G**:PM6) level only. It should be noticed that the most stable geometry depends on the employed functionals, basis sets and semi-empirical methods, as emphasized in the supplementary material. A general trend is that the most stable situation generally corresponds to a A-type site, i.e., above the center of a

ONIOM DFT Calculations
We start by examining the results obtained with the 6 combinations of basis sets, functionals and semi-empirical methods.
As displayed in Table 1, the binding energy of the hydrogen atom appears to be very small, in a consistent pattern for all levels of calculations. The slightly negative values correspond to an artefact due to the convergence threshold of the geometry optimization process. It is possible that using tighter convergence criteria on the geometry optimization process would solve this issue. However, the Potential Energy Surface (PES) is indeed very flat, with the H atom roaming almost freely above the surface (see e.g., Figure 2). Hence the convergence of the SCF process can be very slow, even using the expensive quadratic convergence algorithm (Schlegel and McDouall, 1991) and the computational cost is quite high.
The carbon atom case is also specific since all calculations, with all 6 levels of theory, give a rather large BE, above 1 eV ( Table 1). The length of the C-O bond, is about 1.6 Å (see Figure 3A). This is in agreement with recent theoretical work on crystalline (Ferrero et al., 2020b) and amorphous (Shimonishi et al., 2018) ices who conclude on a chemisorption rather than physisorption. Indeed, in the triplet state of gas phase OCH 2 , the molecule is slightly pyamidalized and the C-O bond is elongated to a larger value of ∼1.88 Å (Hickson et al., 2016). Ferrero et al. also find a transfer of a hydrogen atom to a nearby water molecule to form a COH radical. Furthermore, a Mulliken population analysis shows that this radical has a negative charge. Our own calculations show a slightly different picture: a Wiberg bond index analysis (Wiberg, 1968) confirms a weak covalent bond between C and O (0.73 for C for A2 site). Furthermore, a NBO population analysis (Weinhold and Carpenter, 1988), shows that the COH 2 remains essentially neutral. The only exception is the A4 site, for which, at all levels of calculation, the adsorption of the carbon atom leads to the loss of a H atom from the attacked H 2 O ( Figure 3B). The NBO analysis this time confirms that COH fragment carries a total charge of −0.733 while the H 3 O part has a +0.650 positive charge ( Figure 3C). Thus, this geometry corresponds to the formation of a H 3 O + -COH − complex, very much like in a solvated situation and in agreement with Ferrero et al. (Ferrero et al., 2020b).
On the other hand, the last two studied atoms, namely N and O in their respective ground state, have weak interaction energies, characteristic of physisorption ( Table 1). The most stable geometries are displayed in Figures 4A,B, respectively, Table 1 also shows a remarkably large difference between the M06-2X and ωB97X-D, independent of the basis sets and semiempirical methods used. The M06-2X are about twice larger, despite the fact that the geometries are very similar for all sites. Table 2 summarizes the results for the 4 studied diatomic and triatomic molecules (NH, OH, HCO and H 2 O). The corresponding lowest lying equilibrium geometries are displayed in Figures 5, 6. The maximum BE of NH is about 50% larger than for OH, which is typical of hydrogen bonding (David et al., 2017). For OH, the maximum BE, around 0.7 eV depending on the level of calculation, is typical of the formation of two hydrogen bonds, in agreement with previous theoretical works (Du et al., 2006;Guedes et al., 2003), as shown in Figure 5B.
The situation is very similar in the HCO case, which creates two hydrogen bonds as well ( Figure 6A) although the BE is slightly lower (about 0.5 eV). For H 2 O, the rather high value (about 0.6 eV according to Table 2), is consistent with the building of three hydrogen bonds in the most favorable conformation shown in Figure 6B.
The BE of CH 3 is consistently predicted to be around 0.3 eV ( Table 3) and this radical is only bound by a weak interaction as can be seen in Figure 7A. In contrast, for ammonia, the calculated binding energies appear to be rather high, between 0.8 and 1.2 eV depending on the level of theory (Table 3). This is correlated to the existence of two hydrogen bonds between NH 3 and two water molecules as displayed in Figure 7B. This is consistent with the intermolecular energy obtained for the isolated NH 3 -H 2 O dimer (0.473 eV) (Wakelam et al., 2017). From the inspection of the average energies given in Tables 1-3, and the fact that all methods happened to give very similar geometries, we have selected the ωB97X-D/6-31+G**:PM6 combination to obtain the ZPE correction and apply the single point CBS correction. Indeed, the M06-2X gives too large values for N and O atoms compared to ωB97X-D. This latter functional, being range-separated describes more accurately long-range interactions. Basis sets effects seem to be small and the 6-31+G** one is much less demanding. Finally, PM6 was selected over PM7R8 since this latter method, although available in the Gaussian 16 code and detailed in the PhD work of Throssell (Throssell, 2017), does not seem to have been published yet. Table 4 shows the summary of the results obtained for crystalline ices at our best level of calculations, i.e., CBS/DLPNO-CCSD(T): PM6//ΟΝΙΟΜ(ωB97X-D/6-31+G**:PM6) including scaled ZPE correction calculated at the ΟΝΙΟΜ(ωB97X-D/6-31+G**: PM6) level. A graphical representation is shown in Figure 8. The impact of ZPE on BE's is shown in Figure 8A Except for the case of CH 3 , a linear fit gives a correlation coefficient larger than 0.96. As expected, ZPE effects are larger for polyatomic molecules. For C, N, and O, the BE tends to increase by 3-7 %. For molecules, the ZPE lowers the BE's, the drop varying from 10 to 20% with the notable exception of CH 3 where it reaches 33%. Recently, Ferrero et al. (Ferrero et al., 2020a) have obtained the relation BE0 0.854 BEe using a fragment approach to calculate ZPE effects on crystalline ice. Their set of species does not include any atomic system but contains larger molecules, leading to much larger ZPE's.

CBS Calculations
In addition, single point CBS/DLPNO-CCSD(T) calculations have also some effect on the BE's, as illustrated in Figure 8B. The correction is more erratic, with a bad (< 0.9) correlation coefficient for H, O, CH 3 , and NH 3 . Except for the nitrogen atom where CBS tends to augment the BE by 9%, the other values are lowered, sometimes considerably (up to 55% for HCO for example). More details on ZPE and CCSD(T) corrections for each of the eight sites can be found in the supplementary material file. Table 4 also compares when available (i.e., OH, HCO, H 2 O, and CH 3 ) the present results to the single point DFT calculations on Mechanical Embedding (ME) ONIOM(DFT/AMOEBA) geometries of Sameera et al. (Sameera et al., 2017;Sameera and Maseras, 2018). It should be noticed that since these calculations do not include any ZPE, the BE's are larger than the present results. In addition, the maximum values given in Table 4 correspond to the most stable configuration on a perfect crystalline ice. Our results can also be compared to the recent values obtained by Ferrero et al. (Ferrero et al., 2020a;Ferrero et al., 2020b) using periodic DFT methods and given at the bottom of Table 5. Except for N atom, our values for C and O are much lower than the periodic PBE-D3(BJ)/TZVP of Ferrero et al. (Ferrero et al., 2020b). The recent study of Martinez-Bachs et al. (Ferrero et al., 2020b) performed in an ONIOM-like CCSD(T)// M06-2X/TZVP scheme provides a rather high value for N (0.143 eV) compared to our maximum (0.083 eV). On the other hand, their value for NH (0.364 eV) is lower than our maximum result equal to 0.503 eV. Except for the HCO case where the two values coincide, CBS values are also lower than the CCSD(T) ones obtained by Ferrero et al. (Ferrero et al., 2020a), but it should be kept in mind that these authors limited the calculation to only two water molecules in addition to the adsorbate. Amorphous ice Table 5 shows the summary of the results obtained for amorphous ices at the CBS/DLPNO-CCSD(T):PM6// ΟΝΙΟΜ(ωB97X-D/6-31+G**:PM6) including scaled ZPE correction. This table also displays when available previous theoretical determinations of the BE's. A graphical comparison with the crystalline case can be seen in Figure 11. Figure 9A compares the corrected BE 0 to the uncorrected BE e DFT results. If we except the special case of hydrogen atom, the correlation coefficient is always larger than 0.97. For the three other atoms, the ZPE correction is in fact very small (a few %). For the molecules, the correction is between 11 and 20 %, except for CH 3 where it is 32% as in the crystalline case. Figure 9B plots the comparison between the CBS/DLPNO-CCSD(T) and the DFT calculations. The dispersion of the values is much larger than in the crystalline case, except for the carbon atom where a covalent bond is created. However, it is clear that the CBS results are considerably lower than the DFT ones. It should be noticed that, as in the A4 site of crystalline ice, one of the amorphous clusters produced a H 3 O + COH − complex when approached by a carbon atom in its triplet ground state (see Figure 10A). The charges of this product were confirmed via a NBO analysis ( Figure 10B). The H 3 O + COH − structure lies about 1.1 eV below the dissociated case, which is the largest calculated value. For the seven other amorphous structures, the final product was the pyramidal triplet state of the COH 2 complex already found in crystalline ice. Table 5, our best estimate for the binding energy of H is 0.010 ± 0.020 eV. However, the maximum value is 0.058 eV. They are close to the crystalline case (Table 4). Previous theoretical works, either using quantum methods on clusters [(Wakelam et al., 2017;Das et al., 2018)] or classical calculations (Al-Halabi and Van Dishoeck, 2007) also give very dispersed values. One can only conclude that while on average the BE is very small for H, some local configuration of ASW may lead to larger values as large as 0.058 eV.

According to
Concerning the carbon atom, the low value obtained by Das et al. (Das et al., 2018), i.e., only 0.057 eV is easily explained by the fact that the C atom is weakly bound to a hydrogen atom instead of an oxygen [see Figure 3A in supplementary material in (Das et al., 2018)]. Our average value of 0.816 ± 0.132 eV is 33% lower than average of 1.215 ± 0.036 eV obtained at the ωb97X-D/6-311 + G(d, p) level by Shimonishi et al. (Shimonishi et al., 2018). The differences come from the lack of Low level layer, the absence of ZPE and CBS single point correction which impact the BE's. Their smaller uncertainty is possibly due to a different method to select the adsorption trial geometries and average the largest BE's.
The binding energy of nitrogen atom on ASW has been measured experimentally by Minissale et al. (Minissale et al., 2016) to be 0.062 +0.014 −0.007 eV. This is in reasonable agreement with our CBS result of 0.050 ± 0.020 eV but somewhat larger than the recent values of 0.034 ± 0.003 eV (Shimonishi et al., 2018) and 0.030 ± 0.013 eV (Molpeceres et al., 2020).

Frontiers in Astronomy and
To the best of our knowledge, there is no experimental determination of the BE of the NH radical on ice surfaces. As shown in Table 5, only calculations on water clusters are available. Our best estimate of 0.210 ± 0.097 eV is slightly higher than the 0.168 eV value obtained on a water tetramer by Das et al. (Das et al., 2018) but close to the dimer result (0.207 eV) from Wakelam et al. (Wakelam et al., 2017).
Similarly, for the OH radical, no experimental value has been published yet. For example, for OH,  assumed that it is half the value of H 2 O, i.e., 0.246 eV. Dulieu et al. (Dulieu et al., 2013) obtained 0.396 eV but on silicate dust surfaces. Therefore, as for other molecules, only theoretical values are available, most of them obtained with water clusters ( Table 5). The only exception is the recent work by Ferrero et al. (Ferrero et al., 2020b), where the BE of OH is calculated to be between 0.134 and 0.459 eV (on 5 amorphous structures). Our maximum value of 0.769 eV is higher, and could be explained by the existence of three hydrogen bonds between OH and the neighboring water molecules. Our average best estimate is thus 0.448 ± 0.167 eV.
As for HCO, again no experimental value is available. The value of  of 0.138 eV was deduced by addition and subtraction of energies from other molecules. In Table 5 In contrast to most other species studied in the present work, the adsorption of water on ice has been the subject of numerous experimental works [see Table 1 in (Yocum et al., 2019) for a recent review]. Sandford and Allamandola (Sandford and Allamandola, 1988) proposed 0.415 ± 0.001 and 0.437 ± 0.001 eV, Fraser et al. (Fraser et al., 2001) 0.497 ± 0.005 eV, Haynes et al. (Haynes et al., 1992) 0.516 ± 0.009 eV, Bolina et al. (Bolina et al., 2005) 0.414 ± 0.008 eV, Ulbricht et al. (Ulbricht et al., 2006) 0.477 ± 0.031 eV. Although the experiments differ in the details (for example the substrate), they agree to a value between 0.4-0.5 eV. This is in good agreement with our best estimate of 0.467 ± 0.067 eV ( Table 5).
CH 3 being a radical as well, experimental determinations seem inexistent. A recent DFT study by Enrique-Romero et al.  Table 5). Our own results are slightly higher with 0.114 ± 0.036 eV. Our minimum and maximum values are also similar to those of Ferrero et al. (Ferrero et al., 2020b).
Finally, for NH 3 , the recommended value of 0.477 eV (Hama and Watanabe, 2013) is estimated from the experimental TPD experiments from Collings et al. (Collings et al., 2004). As can be seen in Table 4, this value is close to the water clusters calculation [ (Wakelam et al., 2017;Das et al., 2018)] and lies within the interval of our best estimate: 0.371 ± 0.120 eV. The DFT//HF-3c results of Ferrero et al. (Ferrero et al., 2020b), favor larger values between 0.372 and 0.650 eV. Figure 11. A gives finally a graphical representation of the CBS BEs obtained on the 16 adsorption sites, 8 crystalline and 8 amorphous structures. It shows that BE's can vary significantly Frontiers in Astronomy and Space Sciences | www.frontiersin.org March 2021 | Volume 8 | Article 645243 for a given adsorbate or a given substrate, depending strongly on the local environment. Excepting the reactive carbon atom case, the highest value should correspond to the lowest lying equilibrium geometry. However, in most cases, several configurations have very similar BE's, both for crystalline and amorphous cases, these values being indicative of a configuration that may exist on any ice sample. It is worth mentioning that the statistics could be improved by considering a larger number of initial configurations, but it has to be reminded that for each geometry, a full QM/MM optimization and single point CBS/DLPNO-CCSD(T) correction are carried out, these calculations being expensive, more especially when convergence is difficult to achieve. Figure 11B compares the average BE's and their standard deviation for both surfaces. In most cases, the amorphous average is above the crystalline value, but the trend is not universal and does not allow to conclude without ambiguity. It should be recalled that the considered amorphous ice surface is non-porous, addition of pores or cavities could lead to more marked differences with respect to the crystalline face. However, such a surface would be difficult to handle with the proposed methodology because of the large cluster size it would require. The comparison of our calculated BE's with the values that have been reported in UMIST (McElroy et al., 2013) or KIDA (Wakelam et al., 2012) databases indicate that the values appearing in the databases are for most cases contained in the BE's range we obtain on both surfaces. For OH, HCO, and C, the data implemented in KIDA have a better agreement with our theoretical values than the UMIST ones. This comparison further supports our methodology that can be used to provide BE's for species remaining undocumented, such as radicals or atoms or even larger molecules, thus improving the chemical networks implemented in the astrochemical models.

CONCLUSION
In this work we have used a first principles methodology tailored to study adsorption processes of small atoms and molecules on both crystalline and amorphous ices, represented by a cluster of about 150 water molecules. A full quantum ONIOM QMHigh/QMLow scheme was employed, using semi-empirical methods for the surrounding water molecules representative of the ice sample. Several combinations of functionals, basis sets and semiempirical parametrizations were tested on the different possible adsorption sites of a perfect hexagonal ice. The ΟΝΙΟΜ(ωB97X-D/6-31+G**:PM6) was chosen as a good compromise between cost and accuracy. Despite the fact that the ωB97X-D functional is range-separated and includes dispersion correction, the energy of the QM part was recalculated using single-point CBS/DLPNO-CCSD(T) calculations. This level of theory should remove possible BSSE effects. Inclusion of ZPE effect is also important, even for atoms.
This protocol was applied to both crystalline and amorphous ices. For the latter, eight snapshots obtained from classical MD calculations were selected. This allows us to propose accurate binding energies for the set of studied species. The dispersion of the BE's is rather large, reflecting a large diversity in the hydrogen bond network and the possible inclusion of the adsorbate on the surface. In the case of carbon atom, two different pathways were identified, depending on the local environment. One of them leads to the formation of a covalently bonded COH 2 molecule while the other one, less frequent, produces a COH − H 3 O + complex. The present method is sufficiently versatile to be applied to other systems, including larger adsorbed molecules and/or, mixed ices. The accuracy can also be further benchmarked by adapting the choice of the functional and basis set. A possible extension would be to use a three-layer ONIOM. The computational cost of the DLPNO-CCSD(T) single point calculations can also be limited by using only the def2-TZVPP basis set, for example. It would also be interesting to increase the number of studied snapshots in order to get more statistically significant values. Finally, reactive processes may also be tractable.
As a conclusion, the present study brings new BEs estimations to be included in astrochemical models. When pure quantum chemistry calculations are performed, a single and unique value is produced that corresponds to the global minimum. With the methodology proposed here, the statistics over different configurations provides a set of values: maximum, minimum, average and standard deviation, BE distribution. The average BE accounts for the diversity of adsorption sites, whereas using the maximum BE implies that adsorption occurs preferentially in the most stable sites. An extensive investigation should be carried out to quantify how sensitive the astrochemical models are with respect to BE variations. This is beyond the scope of the present study but would help to target the BE accuracy and consequently to adjust the theoretical methodology accordingly.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
FIGURE 11 | (A) Summary of CBS BE (in eV) distributions for crystalline (blue) and amorphous (orange) ices (B) Comparison of calculated average BE's (with standard deviation) between crystalline (blue) and amorphous (orange) ice. Values, when available, taken from UMIST (McElroy et al., 2013) and KIDA (Wakelam et al., 2012) databases are also represented in grey and yellow, respectively.