<figure><img src="image_1.png" /></figure>

# Binding Energies of Interstellar Molecules on Crystalline and Amorphous Models of Water Ice by Ab Initio Calculations

Stefano Ferrero<sup>1,2 Ⓓ</sup>, Lorenzo Zamirri<sup>2,3 Ⓓ</sup>, Cecilia Ceccarelli<sup>4 Ⓓ</sup>, Arezu Witzel<sup>4 Ⓓ</sup>, Albert Rimola<sup>1 Ⓓ</sup>, and Piero Ugliengo<sup>2 Ⓓ</sup>

Received 2020 April 28; revised 2020 August 25; accepted 2020 September 15; published 2020 November 17

## Abstract

In the denser and colder (≤20 K) regions of the interstellar medium (ISM), near-infrared observations have revealed the presence of submicron-sized dust grains covered by several layers of H₂O-dominated ices and “dirtied” by the presence of other volatile species. Whether a molecule is in the gas or solid-phase depends on its binding energy (BE) on ice surfaces. Thus, BEs are crucial parameters for the astrochemical models that aim to reproduce the observed evolution of the ISM chemistry. In general, BEs can be inferred either from experimental techniques or by theoretical computations. In this work, we present a reliable computational methodology to evaluate the BEs of a large set (21) of astrochemical relevant species. We considered different periodic surface models of both crystalline and amorphous nature to mimic the interstellar water ice mantles. Both models ensure that hydrogen bond cooperativity is fully taken into account at variance with the small ice cluster models. Density functional theory adopting both B3LYP-D3 and M06-2X functionals was used to predict the species/ice structure and their BEs. As expected from the complexity of the ice surfaces, we found that each molecule can experience multiple BE values, which depend on its structure and position at the ice surface. A comparison of our computed data with literature data shows agreement in some cases and (large) differences in others. We discuss some astrophysical implications that show the importance of calculating BEs using more realistic interstellar ice surfaces to have reliable values for inclusion in the astrochemical models.

*Unified Astronomy Thesaurus concepts: Surface ices (2117); Interstellar dust (836); Interstellar molecules (849); Dense interstellar clouds (371); Interstellar medium (847); Solid matter physics (2090); Interstellar dust processes (838); Computational methods (1965)*

*Supporting material: figure set*

## 1. Introduction

The presence of molecules in the extreme physical conditions of the interstellar medium (ISM) was considered impossible by astronomers, until the first diatomic species (CN, CH, and CH⁺) were detected in the ISM from optical and ultraviolet transitions (Swings & Rosenfeld 1937; McKellar 1940; Douglas & Herzberg 1942). Nowadays more than 200 gaseous molecular species (including radicals and ions) have been identified in the diffuse and dense regions of the ISM, thanks to their rotational and vibrational lines in the radio to far-infrared (FIR) wavelengths (e.g., see the review by McGuire 2018). In the coldest (≤20–90 K) and densest (≥10³ cm⁻<fnref n="3" />) ISM, some of these molecules are also detected in the solid state via near-infrared (NIR) observations (e.g., see the review by Boogert et al. 2015).

We now know that the solid-state molecules are frozen species that envelop the submicron dust grains that permeate the ISM and whose refractory core is made of silicates and carbonaceous materials (e.g., Jones 2013; Jones et al. 2017). The grain iced mantle composition is governed by the adsorption of species from the gas phase and by chemical reactions occurring on the grain surfaces. For example, the most abundant component of the grain mantles is H₂O, which is formed by the hydrogenation of O, O₂, and O₃ on the grain surfaces (e.g., Hiraoka et al. 1998; Dulieu et al. 2010; Oba et al. 2012).

The water-rich ice is recognized from two specific NIR bands at about 3 and 6 µm, which are associated with its O–H

stretching and H–O–H bending modes, respectively (e.g., see the review by Boogert et al. 2015). In addition, species like CO, CO₂, NH₃, CH₄, CH₃OH, and H₂CO have also been identified as minor constituents of the ice mantles, which, for this reason, are sometimes referred to as “dirty ices” (Boogert et al. 2015). Furthermore, the comparison between the astronomical spectroscopic observations and the laboratory spectra of an analogous interstellar ice sample, principally based on the O–H stretching feature, has shown that the mantle ices very likely possess an amorphous-like structure resembling that of amorphous solid water (ASW; e.g., Watanabe & Kouchi 2008; Oba et al. 2009; Boogert et al. 2015).

Ice surfaces are known to have an important role in the interstellar chemistry because they can serve as catalysts for chemical reactions that cannot proceed in the gas phase, such as the formation of H₂, the most abundant molecule in the ISM (Hollenbach & Salpeter 1971). Ice surfaces can catalyze reactions by behaving as (i) a passive third body, this way absorbing part of the excess of energy released in the surface processes (adsorption and/or chemical reaction) (e.g., Panta-leone et al. 2020); (ii) a chemical catalyst, this way directly participating in the reaction reducing the activation energies (e.g., Rimola et al. 2018; Enrique-Romero et al. 2019, 2020); or (iii) a reactant concentrator, this way retaining the reactants and keeping them in close proximity for subsequent reaction (e.g., CO adsorption and retention for subsequent hydrogenation to form H₂CO and CH₃OH (e.g., Watanabe & Kouchi 2002; Rimola et al. 2014; Zamirri et al. 2019b). All three processes

<footnote n="1"><sup>1</sup>Departament de Química, Universitat Autonoma de Barcelona, E-08193 Bellaterra, Catalonia, Spain</footnote>

<footnote n="2"><sup>2</sup>Dipartimento di Chimica, Università degli Studi di Torino, via P. Giuria 7, I-10125, Torino, Italy; piero.ugliengo@unito.it</footnote>

<footnote n="3"><sup>3</sup>Nanostructured Interfaces and Surfaces (NIS) Centre, Università degli Studi di Torino, via P. Giuria 7, I-10125, Torino, Italy</footnote>

<footnote n="4"><sup>4</sup>Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France; cecilia.ceccarelli@univ-grenoble-alpes.fr</footnote>
depend on the binding energies (BEs) of the molecules either directly (e.g., the adsorption of the species) or indirectly (e.g., because the diffusion of a particle on the grain surfaces is a fraction of its BE) (see Cuppen et al. 2017). In addition, molecules formed on the grain surfaces can be later transferred to the gas phase by various desorption processes, most of which depend, again, on the BE of the species. In practice, BEs are crucial properties of the interstellar molecules and play a huge role in the resulting ISM chemical composition. This key role of BEs is very obvious in the astrochemical models that aim at reproducing the chemical evolution of interstellar objects, as clearly shown by two recent works by Wakelam et al. (2017) and Penteado et al. (2017), respectively.

Experimentally, the BEs of astrochemical species are measured by temperature programmed desorption (TPD) experiments. These experiments measure the energy required to desorb a particular species from the substrate, namely, a desorption enthalpy, which is equal to the BE only if there are no activated processes (He et al. 2016) and if thermal effects are neglected. A typical TPD experiment consists of two phases. In the first one, the substrate, maintained at a constant temperature, is exposed to the species that have to be adsorbed coming from the gas phase. In the second phase, the temperature is increased until desorption of the adsorbed species—collected and analyzed by a mass spectrometer—occurs. The BE is then usually extracted by applying the direct inversion method on the Polanyi–Wigner equation (e.g., Dohnalek et al. 2001; Noble et al. 2012). The BE values obtained in this way strongly depend on the chemical composition and morphology of the substrate and also on whether the experiment is conducted in the monolayer or multilayer regime (e.g., Noble et al. 2012; He et al. 2016; Chaabouni et al. 2018). Another issue related to the TPD technique is that it cannot provide accurate BBEs for radical species, as they are very reactive. In the literature, there are many works that have investigated the desorption processes by means of the TPD technique (e.g., Collings et al. 2004; Noble et al. 2012; Dulieu et al. 2013; Fayolle et al. 2016; He et al. 2016; Smith et al. 2016), but they have been conducted for just a handful of important astrochemical species, whereas a typical network of an astrochemical model can contain up to 500 species and very different substrates. In a recent work, Penteado et al. (2017) collected the results of these experimental works, trying to be as homogeneous as possible in terms of different substrates, estimating the missing BE values from the available data and performing a systematic analysis on the effect that the BE uncertainties can have on astrochemical model simulations.

BE values can also be obtained by means of computational approaches that, in some situations, can overcome the experimental limitations. Many computational works have so far focused on a few important astrochemical species like H, H₂, N, O, CO, and CO₂, in which BEs are calculated on periodic/cluster models of crystalline/amorphous structural states using different computational techniques (e.g., Al-Halabi & Van Dishoeck 2007; Karssemeijer et al. 2014; Karssemeijer & Cuppen 2014; Ásgeirsson et al. 2017; Senevirathne et al. 2017; Shimonishi et al. 2018; Zamirri et al. 2019a). In addition, other works have computed BEs in a larger number of species but with a very approximate model of the substrate. For example, in a recent work by Wakelam et al. (2017) BE values of more than 100 species are calculated by approximating the ASW surface with a single water molecule. The authors then fitted the most reliable BE measurements (16 cases) against the corresponding computed ones, obtaining a good correlation

between the two data sets. In this way, all the errors in the computational methods and limitations due to the adoption of a single water molecule are compensated by the fitting with the experimental values, in the view of the authors. The resulting parameters are then used to scale all the remaining computed BEs to improve their accuracy. This clever procedure does, however, consider the proposed scaling universal, leaving aside the complexity of the real ice surface and the specific features of the various adsorbates. In a similar work, Das et al. (2018) have calculated the BEs of 100 species by increasing the size of a water cluster from one to six molecules, noticing that the calculated BE approaches the experimental value when the cluster size is increased. As we will show in the present work, these approaches, relying on an arbitrary and very limited number of water molecules, cannot, however, mimic a surface of icy grain. Furthermore, the strength of interaction between icy water molecules, as well as with respect to the adsorbates, depends on the hydrogen bond cooperativity, which is underestimated in small water clusters.

In this work, we followed a different approach, focusing on extended periodic ice models, either crystalline or amorphous, adopting a robust computational methodology based on a quantum mechanical approach. We simulate the adsorption of a set of 21 interstellar molecules, 4 of which are radical species, on several specific exposed sites of the water surfaces of both extended models. BE values have been calculated for more than one binding site (if present) to provide the spread of the BE values that the same molecule can have depending on the position in the ice. Different approaches, with different computational cost, have been tested and compared, and the final computed BEs have been compared with data from the computational approaches of Wakelam et al. (2017) and Das et al. (2018) and data from UMIST and KIDA databases, as well as available experimental data (e.g., McElroy et al. 2013; Wakelam et al. 2015). One added value of this work is the definition of both a reliable, computationally cost-effective ab initio procedure designed to arrive at accurate BE values and an ice grain atomistic model that can be applied to predict the BEs of any species of astrochemical interest.

## 2. Computational Details

### 2.1. *Structure of the Ice: Periodic Simulations*

Water ice surfaces have been modeled enforcing periodic boundary conditions to define icy slabs of finite thickness either entirely crystalline or of amorphous nature. Adsorption is then carried out from the void region above the defined slabs. Periodic calculations have been performed with the ab initio CRYSTAL17 code (Dovesi et al. 2018). This software implements both the Hartree–Fock (HF) and Kohn–Sham self-consistent field methods for the solution of the electronic Schrödinger equation, fully exploiting, if present, the crystal–line or molecular symmetry of the system under investigation. CRYSTAL17 adopts localized Gaussian functions as basis sets, similar to the approach followed by molecular codes. This allows CRYSTAL17 to perform geometry optimizations and vibrational properties of both periodic (polymer, surfaces, and crystals) and nonperiodic (molecules) systems with the same level of accuracy. Furthermore, the definition of the surfaces through the slab model allows us to avoid the 3D fake replica of the slab as forced when adopting the plane wave basis set.
Computational parameters are set to values ensuring good accuracy in the results. The threshold parameters for the evaluation of the Coulomb and exchange bi-electronic integrals (TOLINTEG keyword in the CRYSTAL17 code; Dovesi et al. 2018) have been set equal to 7, 7, 7, 7, and 14. The needed density functional integration is carried out numerically over a grid of points, which is based on an atomic partition method developed by Becke (1988). The standard pruned grid (XLGRID keyword in the CRYSTAL17 code; Dovesi et al. 2018), composed of 75 radial points and a maximum of 974 angular points, was used. The sampling of the reciprocal space was conducted with a Pack–Monkhorst mesh (Pack & Monkhorst 1977), with a shrinking factor (SHRINK in the code CRYSTAL17; Dovesi et al. 2018) of 2, which generates 4k points in the first Brillouin zone. The choice of the numerical values we assigned to these three computational parameters is fully justified in Appendix A.

Geometry optimizations have been carried out using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm (Broyden 1970; Fletcher 1970; Goldfarb 1970; Shanno 1970), relaxing both the atomic positions and the cell parameters. We adopted the default values for the parameters controlling the convergence, i.e., difference in energy between two subsequent steps, $1 \times 10^{-7}$ Hartree; and maximum components and rms of the components of the gradients and atomic displacement vectors, $4.5 \times 10^{-4}$ hartrees bohr$^{-1}$ and $3 \times 10^{-4}$ hartrees bohr$^{-1}$, and $1.8 \times 10^{-3}$ bohr and $1.2 \times 10^{-3}$ bohr, respectively. All periodic calculations were grounded on either the density functional theory (DFT) or the HF-3c method (Hohenberg & Kohn 1964; Sure & Grimme 2013). Within the DFT framework, different functionals were used to describe closed- and open-shell systems. For the former, we used the hybrid B3LYP method (Lee et al. 1988; Becke 1993), which has been shown to provide a good level of accuracy for the interaction energies of noncovalent bound dimers (Kraus & Frank 2018), added with the D3-BJ correction for the description of dispersive interactions (Grimme et al. 2010, 2011). For open-shell systems, treated with a spin-unrestricted formalism (Pople et al. 1995), we used the hybrid M06-2X functional (Zhao & Truhlar 2008), which has been proved to give accurate results in estimating the interaction energy of noncovalent binary complexes involving a radical species and a polar molecule (Tentscher & Arey 2013). The choice of these two different functionals is justified by two previous works describing the accuracy on the energetic properties of molecular adducts (Tentscher & Arey 2013; Kraus & Frank 2018). For all periodic DFT calculations we used the Ahlrichs triple-zeta quality VTZ basis set, supplemented with a double set of polarization functions (Schäfer et al. 1992). In the following, we will refer to this basis set as “A-VTZ” (see Appendix D for details of the adopted basis set).

The HF-3c method is a new method combining the Hartree–Fock Hamiltonian with the minimal basis set MINI-1 (Tatewaki & Huzinaga 1980) and with three a posteriori corrections for (i) the basis set superposition error (BSSE), arising when localized Gaussian functions are used to expand the basis set (Jansen & Ros 1969; Liu & McLean 1973); (ii) the dispersive interactions; and (iii) short-ranged deficiencies due to the adopted minimal basis set (Sure & Grimme 2013).

Harmonic frequency calculations were carried out on the optimized geometries of both crystalline and amorphous ices to characterize the stationary points of each structure. Vibrational frequencies have been calculated at the Γ point by diagonalizing

the mass-weighted Hessian matrix of second-order energy derivatives with respect to atomic displacements (Pascale et al. 2004; Zicovich-Wilson et al. 2004). The Hessian matrix elements have been evaluated numerically by a six-point formula (NUMDERIV = 2 in the CRYSTAL17 code; Dovesi et al. 2018), based on two displacements of ±0.001 Å for each nuclear Cartesian coordinate from the minimum structure.

To avoid computational burden, only a portion of the systems has been considered in the construction of the Hessian matrix, including the adsorbed species and the spatially closest interacting water molecules of the ice surface. This “fragment” strategy for the frequency calculation has already been tested by some of us in previous works and is fully justified by the noncovalent nature of the interacting systems, where the coupling between the vibrational modes of bulk ice and adsorbate moieties is negligible (Tosoni et al. 2005; Rimola et al. 2008; Zamirri et al. 2017).

From the set of frequencies resulting from the “fragment” calculations we worked out the zero-point energy (ZPE) for the free crystalline ice surface, the free adsorbate, and the ice surface/adsorbate complex to arrive at the corresponding correction ΔZPE, as reported in Appendix A.1. From the ΔZPE we corrected the electronic BE for each adsorbate as $BE(0) = BE - \Delta ZPE$ and found a good linear correlation $BE(0) = 0.854 BE$, as shown in Appendix A.1. While the “fragment frequency” strategy is fine for computing the ΔZPE of the crystalline ice model owing to the structural rigidity enforced by the system symmetry, the same does not hold for the amorphous ice. In that case, the large unit cell (60 water molecules) and their random organization render the ice structure rather sensitive to the adsorbate interaction, which causes large structural water molecule rearrangement. This, in turn, alters significantly the whole set of normal modes, and the numerical value of the ΔZPE becomes ill-defined. Nevertheless, considering that the kinds of interactions operative for the crystal ice are of the same nature as those for the amorphous one, we adopted the same scaling factor of 0.854 computed for the crystalline ice to correct the electronic BE for the amorphous one. In the following, we compared the experimental BE usually measured for amorphous ices with the $BE(0)$ values. To discuss the internal comparison between adsorption features of different adsorbates on the crystalline ice, we still focused on the uncorrected BEs.

## 2.2. BE Calculation and Counterpoise Correction

When Gaussian basis sets are used, a spurious contribution arises in the calculation of the molecule/surface interactions, called the BSSE (e.g., Boys & Bernardi 1970). In this work, the BSSE for DFT calculations has been corrected making use of the a posteriori counterpoise (CP) correction by Boys and Bernardi (Davidson & Feller 1986). The CP-corrected interaction $\Delta E^{CP}$ energy has been calculated as

$$ \Delta E^{CP} = \Delta E^* + \delta E + \Delta E_L - \text{BSSE} \quad (1) $$

where $\Delta E^*$ is the deformation-free interaction energy, $\delta E$ is the total contribution to the deformation energy, and $\Delta E_L$ is the lateral interaction (adsorbate–adsorbate interaction) energy contribution. Details on the calculation of each energetic term of Equation (1) can be found in Appendix A. By definition, BE
is the opposite of the CP-corrected interaction energy:

$$ \Delta E^{\text{CP}} = -BE. \qquad (2) $$

## 2.3. BE Refinement with the Embedded Cluster Method

With the aim of refining the periodic DFT BE values for the crystalline ice model, single-point energy calculations have been carried out on small clusters, cut out from the crystalline ice model, using a higher level of theory than the DFT methods with the Gaussian09 program (Frisch et al. 2009). The adopted cluster models were derived from the periodic systems and are described in Section 3.2.2. These refinements have been performed through the ONIOM2 approach (Dapprich et al. 1999), dividing the systems into two parts that are described by two different levels of theory. The Model system (i.e., a small moiety of the whole system, including the adsorbate and the closest water molecules) was described by the High level of theory represented by the single- and double-electronic excitation coupled-cluster method with an added perturbative description of triple excitations (CCSD(T)). The Real system (i.e., the whole system) was described by the DFT level of theory adopted in the periodic calculations with the two different functionals for open- and closed-shell species. In the ONIOM2 methodology, the BE can be written as

$$ \text{BE}(\text{ONIOM2}) = -\text{BE}(\text{Low, Real}) + \Delta\text{BE} \qquad (3) $$

$$ \Delta\text{BE} = \text{BE}(\text{High, Model}) - \text{BE}(\text{Low, Model}). \qquad (4) $$

The final BE(ONIOM2) is also corrected by the BSSE following the same scheme described above. Our choices about the Model and Real systems will be extensively justified in Section 3.2.2.

# 3. Results

## 3.1. Ice Surface Models

### 3.1.1. Crystalline Ice Model

Despite the amorphous and perhaps porous nature of the interstellar ice, we adopted, as a paradigmatic case, a proton-ordered crystalline bulk ice model usually known as P-ice (Pna2₁ space group; Casassa et al. 1997). From P-ice bulk, we cut out a slab model, i.e., a 2D-periodic model representing a surface. Consequently, periodic boundary conditions are maintained only along the two directions defining the slab plane, while the third direction (z-axis) is nonperiodic and defines the slab thickness. The slab model adopted in this work represents the P-ice (010) surface, in accordance with previous work (Zamirri et al. 2018). This slab consists of 12 atomic layers, is stoichiometric, and has a null electric dipole moment across the z-axis. This ensures an electronic stability of the model with the increase of the slab thickness (Tasker 1979). The slab structure has been fully optimized (unit cell and atomic fractional coordinates) at both B3LYP-D3/A-VTZ* and M06-2X/A-VTZ* DFT levels. As can be seen from Figure 1 (panel (a)), the (010) P-ice unit cell is rather small, showing only one dangling hydrogen (dH) and oxygen (dO) as binding sites. For large molecules, to increase the number of adsorption sites and minimize the lateral interactions among replicas of the adsorbate, we also considered a 2 × 1 supercell. The electrostatic potential maps (EPMs; see Figure 1, panels (b) and (c)) clearly reveal positive (blue EPM regions) and

<figure><img src="image_2.png" /><figcaption>Figure 1. The (010) slab model of P-ice. (a) Side view along the *b* lattice vector. (b) Top view of the 1 × 1 unit cell (|*a*| = 4.500 Å and |*b*| = 7.078 Å) superimposed onto the EPM. (c) Top view of the 2 × 1 supercell (|*a*| = 8.980 Å, |*b*| = 7.081 Å) superimposed onto the EPM. The isosurface value for the electron density where the electrostatic potential is mapped is set equal to 10⁻⁶ au. Color code: +0.02 au (blue, positive), 0.00 au (green, neutral), and −0.02 au (red, negative).</figcaption></figure>

negative (red EPM regions) potentials around the dH and dO sites, respectively.

### 3.1.2. Amorphous Solid Water (ASW)

As anticipated, the (010) P-ice surface might not be a physically sound model to represent actual interstellar ice surfaces, due to the evidence, from the spectroscopic feature of the interstellar ice, of its amorphous nature (Boogert et al. 2015). The building up of amorphous surface models is a nontrivial and not unique procedure, because of the lack of a consistent and universally accepted strategy. One common approach is to start from a crystalline model and heat it up to relatively high temperature by running molecular dynamics (MD) simulations for a few picoseconds. This step is followed by thermal annealing to freeze the ice in a glassy amorphous state. In this work, we adopted a different strategy. We refer to a recent work by Shimonishi et al. (2018) in which the BEs of a set of atomic species were computed on several water clusters, previously annealed with MD simulations. We reoptimized (at the B3LYP-D3/A-VTZ* level only) the whole set of ice clusters, and the three most stable clusters, composed of 20 water molecules each, were merged together to define a unit cell of an amorphous periodic ice. This procedure mimics somehow the collision of nanometric-scale icy grains occurring in the molecular clouds. The merger of the three clusters was carried out by matching the dH regions of one cluster with the dO ones of the other. As a result, we ended up with a large 3D-periodic unit cell (with lattice parameters |*a*| = 21.11 Å, |*b*| = 11.8 Å, and |*c*| = 11.6 Å) envisaging 60 water molecules. This initial bulk model was optimized at HF-3c level in order to fully relax the structure from the internal tensions of the initial guess. After this step, we cut out a 2D-periodic slab from the bulk structure. The amorphous slab is composed of 60 water molecules in the unit cell and was further fully optimized (unit cell size and atomic coordinates) at the HF-3c, B3LYP-D3/A-VTZ*, and M06-2X/A-VTZ* levels of theory. The three final structures show little differences in the positions of specific
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Side view of the amorphous slab models. The cell parameter *a* is highlighted as a blue line. Electric dipole moments *μ*<sub>e</sub> along the *z*-direction are shown on the right side.</figcaption></figure>

water molecules, and, on the whole, the structures are very similar (Figure 2). The computed electric dipole moment across the nonperiodic direction (1, 0, 0) for the HF-3c, B3LYP-D3, and M06-2X structures, respectively, showed a very good agreement between different models, also considering the dependence of the dipole value on the adopted quantum mechanical method. These amorphous slab models show different structural features for the upper and lower surfaces, which imparts the residual dipole moment across the slab, and consequently exhibit a variety of different binding sites for adsorbates. To characterize the electrostatic features of these sites, which in turn dictate the adsorption process, we resorted to the EPMs for the top/bottom surfaces of each optimized slab (Figure 3). The general characteristics are very similar for the three models, with B3LYP-D3 and M06-2X giving the closest maps. HF-3c tends to enhance the differences between positive and negative regions owing to overpolarization of the electron density caused by the minimal basis set. “Top” surfaces show a hydrophobic cavity (the central greenish region, Figure 3, absent in the P-ice slab, surrounded by dH positive spots). “Bottom” surfaces show several prominent negative regions (from five dOs) mixed with less prominent positive potentials (due to four buried dHs).

## 3.2. BEs on Crystalline Ice

### 3.2.1. BE Computed with DFT/DFT Method

In this work, we simulated the adsorption of 17 closed-shell species and 4 radicals, shown in Figure 4. For each molecule/surface complex, geometry optimizations (unit cell plus all atomic coordinates without constraints) were performed. Initial

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Color-coded EPMs mapped to the electron density for the “top” and “bottom” surfaces of the amorphous slab (HF-3c, B3LYP-D3, and M06-2X optimized geometries). The dO and dH sites are also labeled. The isosurface value for the electron density is set equal to $10^{-6}$ au, to which the electrostatic potential is mapped out. EPM color code: +0.02 au (blue, positive), 0.00 au (green, neutral), and -0.02 au (red, negative).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Set of molecular and radical species adopted within this work for the calculation of BE on different ice models. $O_2$ is an open-shell (spin-triplet) species (Borden et al. 2017).</figcaption></figure>

structures were guessed by manually setting the maximum number of H-bonds between the two partners. The pure role of dispersion is estimated by extracting the D3 contribution from the total energy at the B3LYP-D3 level of theory. The energetics of the adsorption processes were then computed according to Equation (1).

As can be seen from the results of Table 1, a range of interactions of different strength is established between the adsorbed species and the crystalline P-ice surface. Some molecules do not possess a net electric dipole moment, while exhibiting relevant electric quadrupole moments (i.e., $H_2$, $N_2$, and $O_2$) or multipole moments of higher order (i.e., $CH_4$; see their EPMs in Appendix B). For these cases, only weak interactions are established so that BEs are lower than 1800 K (see BE disp values in Table 1). Interestingly, for the $N_2$, $O_2$, and $CH_4$ cases, interactions are almost repulsive if dispersive contributions are not accounted for in the total BE (compare
<table><caption>**Table 1**
Summary of the BE Values (in Kelvin) Obtained for the Crystalline P-ice (010) Slab with DFT//DFT and DFT//HF-3c Methods</caption><thead><tr><th rowspan="2">Species</th><th colspan="3">(010) P-ice Crystalline Slab DFT//DFT</th><th colspan="3">(010) P-ice Crystalline Slab DFT//HF-3c</th></tr><tr><th>BE disp</th><th>BE no disp</th><th>-disp(%)</th><th>BE disp</th><th>BE no disp</th><th>-disp(%)</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>1191</td><td>565</td><td>625(53)</td><td>926</td><td>241</td><td>686(74)</td></tr><tr><td>O<sub>2</sub></td><td>1022</td><td>-373</td><td>1034(137)</td><td>794</td><td>-84</td><td>878(110)</td></tr><tr><td>N<sub>2</sub></td><td>1564</td><td>-72</td><td>1636(104)</td><td>1455</td><td>-180</td><td>1636(160)</td></tr><tr><td>CH<sub>4</sub></td><td>1684</td><td>-229</td><td>1912(113)</td><td>1912</td><td>-349</td><td>2261(118)</td></tr><tr><td>CO</td><td>2357</td><td>698</td><td>1660(71)</td><td>1948</td><td>60</td><td>1888 (97)</td></tr><tr><td>CO<sub>2</sub></td><td>3440</td><td>1540</td><td>1900(55)</td><td>3007</td><td>938</td><td>2069(69)</td></tr><tr><td>OCS</td><td>3476</td><td>120</td><td>3356(97)</td><td>3187</td><td>265</td><td>2923(92)</td></tr><tr><td>HCl</td><td>6507</td><td>4402</td><td>2093(32)</td><td>6314</td><td>3488</td><td>2237(39)</td></tr><tr><td>HCN</td><td>5124</td><td>3067</td><td>2057(29)</td><td>5725</td><td>3271</td><td>3043(48)</td></tr><tr><td>H<sub>2</sub>O</td><td>8431</td><td>6844</td><td>1588(19)</td><td>8431</td><td>6808</td><td>1612(19)</td></tr><tr><td>H<sub>2</sub>S</td><td>5677</td><td>3380</td><td>2297(40)</td><td>5232</td><td>3199</td><td>2105(40)</td></tr><tr><td>NH<sub>3</sub></td><td>7373</td><td>5533</td><td>1852(25)</td><td>7301</td><td>5484</td><td>1816(25)</td></tr><tr><td>CH<sub>3</sub>CN</td><td>7553</td><td>4450</td><td>3103(41)</td><td>6916</td><td>3259</td><td>2598(44)</td></tr><tr><td>CH<sub>3</sub>OH</td><td>8684</td><td>6014</td><td>2670(31)</td><td>8648</td><td>6026</td><td>2237(27)</td></tr><tr><td>H<sub>2</sub>CO-SC1</td><td>5869</td><td>3885</td><td>1985(34)</td><td>5773</td><td>4053</td><td>2369(37)</td></tr><tr><td>H<sub>2</sub>CO-SC2</td><td>6375</td><td>3692</td><td>2682(42)</td><td>6423</td><td>3716</td><td>2057(36)</td></tr><tr><td>HCONH<sub>2</sub>-SC1</td><td>9610</td><td>6459</td><td>3151(33)</td><td>9321</td><td>6158</td><td>3163(34)</td></tr><tr><td>HCONH<sub>2</sub>-SC2</td><td>10079</td><td>6483</td><td>3608(36)</td><td>9634</td><td>6074</td><td>3560(37)</td></tr><tr><td>HCOOH</td><td>9526</td><td>7325</td><td>2189(23)</td><td>9297</td><td>7168</td><td>2117(23)</td></tr><tr><td>HCOOH-SC</td><td>9442</td><td>7301</td><td>2021(21)</td><td>9405</td><td>7541</td><td>1864(20)</td></tr><tr><td>OH•</td><td>6543<sup>a</sup></td><td></td><td></td><td>6795<sup>a</sup></td><td></td><td></td></tr><tr><td>HCO•</td><td>3476<sup>a</sup></td><td></td><td></td><td>3548<sup>a</sup></td><td></td><td></td></tr><tr><td>CH<sub>3</sub>•</td><td>2562<sup>a</sup></td><td></td><td></td><td>2598<sup>a</sup></td><td></td><td></td></tr><tr><td>NH<sub>2</sub>•</td><td>6038<sup>a</sup></td><td></td><td></td><td>6050<sup>a</sup></td><td></td><td></td></tr></tbody></table>

**Notes.** Legend: “BE disp” = BE value including the D3 contribution; “BE no disp” = BE values not including D3 contribution; “-disp(%)” = absolute (percentage) contribution of dispersive forces to the total BE disp.

BE disp with BE no disp values of Table 1). Therefore, the adsorption is dictated by dispersive forces, which counter-balance the repulsive electrostatic interactions. For the H₂ case, electrostatic interactions are attractive mainly because of the synergic effect of both the surface dH and the dO on the negative and positive parts of the H₂ quadrupole, respectively (see Appendix B).

CO, OCS, and CO₂ also exhibit a quadrupole moment, but due to the presence of heteroatoms in the structure, they can also establish H-bonds with the dH site. Consequently, BEs are larger than the previous set of molecules (i.e., >2400 K; see Table 1). For these three cases, pure electrostatic interactions are attractive, but the dispersion contribution is the most dominant one over the total BE values (compare BE disp with BE no disp values of Table 1). CO, in addition to a net quadrupole, also possesses a weak electric dipole, with the negative end at the carbon atom (see its EPM in Figure 16; see also Zamirri et al. 2017). Thus, although the two negative poles (C and O atoms) of the quadrupole can both interact with the positive dH site, the interaction involving the C atom is energetically slightly favored over the O atom (Zamirri et al. 2017, 2019a). Accordingly, we only considered the C-down case, the computed BE being in good agreement with previous works (Zamirri et al. 2017, 2018). OCS also possesses a dipole and can interact with the surface through either its S- or O-ends, through dO or dH sites. However, due to the softer basic character of S compared to O, the interaction through oxygen is preferred and only considered here.

NH₃, H₂O, HCl, HCN, and H₂S are all amphiprotic molecules that can serve as both acceptors and donors of H-bonds from/to

the dH and dO sites. The relative strong H-bonds with the surface result in total BE values that are almost twice as high as the values of the previous set of molecules (i.e., CO, OCS, and COS). Although also in these cases dispersive forces play an important contribution to the BE, the dominant role is dictated by the H-bonding contribution.

For the adsorption of CH₃OH, CH₃CN, and the three carbonyl-containing compounds, i.e., H₂CO, HCONH₂, and HCOOH, all characterized by large molecular sizes, we adopted the 2 × 1 supercell (shown in Figure 1) to minimize the lateral interactions between adsorbates. Consequently, two dHs and two dOs are available for adsorption. Therefore, for some of these species (i.e., the carbonyl-containing ones), we started from more than one initial geometry to improve a better sampling of the adsorption features on the (010) P-ice surface (the different cases on the supercell are labeled as SC1 and SC2 in Table 1, and the geometries are reported in Figure 16). The BE values of these species are among the highest ones, due to the formation of multiple H-bonds with the slab (and therefore increasing the electrostatic contribution to the interactions) and a large dispersion contribution due to the larger sizes of these molecules with respect to the other species.

The adsorption study has also been extended to four radicals (i.e., OH• NH₂• CH₃• HCO•), since they are of high interest owing to their role in the formation of interstellar compounds (Sorrell 2001; Bennett & Kaiser 2007). OH• and NH₂• form strong H-bonds with the dH and dO sites of the slab, at variance with CH₃• and HCO• cases, as shown by the higher BE values. Because of the nature of the M06-2X functional, we cannot separate the dispersion contributions to the total BEs.

<footnote>a For radical species (energy at M06-2X level) we cannot discern between disp and no disp data.</footnote>
Interestingly, in all cases we did not detect transfer of the electron spin density from the radicals to the ice surface, i.e., the unpaired electron remains localized on the radical species upon adsorption.

### 3.2.2. The ONIOM2 Correction and the Accuracy of the DFT//DFT BE Values

As described in Appendix A, the ONIOM2 methodology has been employed to check the accuracy of the B3LYP-D3/A-VTZ\* and M06-2X/A-VTZ\* theory levels, both representing the Low level of calculation. For this specific case, to reduce the computational burden, we only considered 15 species, leaving aside N₂, O₂, H₂O, CH₄, CH₃CN, and CH₃\* radical. Here, the Real system is the periodic P-ice slab model without adsorbed species. Therefore, the BE(Low,Real) term in Equation (3) corresponds to the BEs at the DFT theory levels, hereafter referred to as BE(DFT, Ice). The Model system is carved from the optimized geometry of the periodic system: it is composed of the adsorbed molecule plus *n* (*n* = 2, 6; the latter only for the H₂ case) closest water molecules of the ice surface to the adsorbates. For the Model systems, two single-point energy calculations have been carried out: one at the High level of theory, i.e., CCSD(T), calculated with Gaussian09, and the other at the Low level of theory, employing the same DFT methods as in the periodic calculations, calculated with CRYSTAL17. For the sake of clarity, we renamed the two terms BE(High,Model) and BE(Low,Model) in Equation (3) for any molecular species *μ* as BE(CCSD(T), *μ*-nH₂O) and BE(DFT, *μ*-nH₂O), respectively.

As CCSD(T) is a wave-function-based method, the associated energy strongly depend on the quality of the adopted basis set (Cramer 2002). Consequently, accurate results are achieved only when complete basis set extrapolation is carried out (Cramer 2002); accordingly, we adopted correlation-consistent basis sets (Dunning 1989), here named as cc-pVNZ, where “cc” stands for correlation consistent and *N* stands for double (D), triple (T), quadruple (Q), etc. Therefore, we performed different calculations improving the quality of the basis set from Jun-cc-pVDZ to Jun-cc-pVQZ (and even Jun-cc-pV5Z when feasible) (Bartlett & Musiał, 2007; Papajak et al. 2011), extrapolating the BE(CCSD(T), *μ*-nH₂O) values for *N* → ∞. Figure 5 shows, using NH₃ as an illustrative example, the plot of the BE(CCSD(T), *μ*-nH₂O) values as a function of 1/*L*³, where *L* is the cardinal number corresponding to the *N* value for each correlation-consistent basis set. For all other species, we observed similar trends. This procedure was used in the past to extrapolate the BE value of CO adsorbed at the Mg(001) surface (Ugliengo & Damin 2002).

The procedure gives for the extrapolated BE(CCSD(T), *μ*-nH₂O) a value of 4089 K, in excellent agreement with the value computed by the plain B3LYP-D3/A-VTZ\* at periodic level of 4390 K (see Figure 5). Very similar agreement was computed for all considered species as shown in Figure 6, in which a very good linear correlation is seen between BE(ONIOM2) and BE(DFT). Therefore, we can confidently assume that the periodic B3LYP-D3/A-VTZ\* (closed-shell molecules) or the M06-2X/A-VTZ\* (radical species) plain BE values are reliable and accurate enough and are those actually used in this work.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** BE(X, *μ*-nH₂O) extrapolated value at infinite basis set for the case of NH₃. The dashed-dotted blue line represents the BE computed for the BE(DFT, *μ*-nH₂O) at the DFT//A-VTZ\* level (4390 K). The solid red line represents the linear fit of the BE(CCSD(T), *μ*-nH₂O) values (red squares) calculated with DZ, TZ, QZ, and 5Z basis sets. The extrapolated BE(X, *μ*-nH₂O) at infinite basis set is highlighted in red in the fitting equation (4089 K).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Linear fit between periodic DFT/A-VTZ\* BE values (BE(DFT)) and the basis set extrapolated ONIOM2 BE values (BE(ONIOM2)). All values are in K. Fit parameters are also reported. Legend: 1—H₂; 2—CO; 3—CO₂; 4—HCO\*; 5—OCS; 6—H₂S; 7—HCN; 8—NH₂\*; 9—H₂CO; 10—HCl; 11—OH\*; 12—NH₃; 13—CH₃OH; 14—HCOOH; 15—HCONH₂.</figcaption></figure>

### 3.2.3. BE Computed with Composite DFT/HF-3c Method

In the previous section we proved the DFT/A-VTZ\* as a reliable and accurate method to compute the BEs of molecules and radicals on the crystalline (010) P-ice ice slab. However, this approach can become very computationally costly when moving from the crystalline to amorphous model of the interstellar ice, as larger unit cells are needed to enforce the needed randomness in the water structure. Therefore, we tested the efficiency and accuracy of the cost-effective computational HF-3c method (see Section 2).

To this end, we adopted a composite procedure that has been recently assessed and extensively tested in the previous work by some of us on the structural and energetic features of molecular crystals, zeolites, and biomolecules (Cutini et al. 2016, 2017, 2019). We started from the DFT/A-VTZ\* optimized structure just
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Linear fit between the BE values calculated with the full DFT computational scheme and the BE values calculated with the composite DFT//HF-3c computational scheme for the crystalline ice model (all values in K). Black filled and open circles stand for open-shell and closed-shell species, respectively.</figcaption></figure>

discussed for the crystalline ice. We reoptimize each structure at the HF-3c level to check the changes in the structures resulting from the more approximated method. Then, we run a single-point energy calculation at the DFT/A-VTZ\* (B3LYP-D3 and M06-2X) levels to evaluate the final BE values. The results obtained are summarized in Figure 7, showing a very good linear correlation between the BE values computed as described.

The largest percentage differences are found for the smallest BEs, that is, those dominated by dispersion interactions or very weak quadrupolar interactions (i.e., N₂, O₂, H₂, and CH₄) in which the deficiencies of the minimal basis set encoded in the HF-3c cannot be entirely recovered by the internal corrections. For higher BE values, the match significantly improves, in some cases being almost perfect. Even for radicals, the composite approach gives good results. It is worth mentioning that HF-3c optimized geometries are very similar to the DFT-optimized ones (only slight geometry alterations occurred), indicating that the adducts are well-defined minima in both potential energy surfaces. This successful procedure calibrated on crystalline ice is therefore adopted to model the adsorption of all 21 species on the proposed amorphous slab model, a task that would have been very expensive at the full DFT/A-VTZ\* level.

### 3.3. BEs on Amorphous Solid Water

On the ASW model, due to the presence of different binding sites, a single BE value is not representative of the whole adsorption process as is the case for almost all adsorbates on the crystalline surface. Therefore, we computed the BE with the composite DFT//HF-3c procedure (see Section 3.2.3) by sampling different adsorption sites at both the “top” and “bottom” surfaces of the amorphous slab. The starting initial structures of each adsorbate were set up by hand, following the maximum electrostatic complementarity between the EPMs (see Figure 3) of the ice surface and that of a given adsorbate. For each molecule, at least four BE values have been computed on different surface sites. Figure 8 reports the examples of methanol and formamide: for each molecule, we show the geometry on the crystalline ice and in two different sites of

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of the final optimized geometries for CH₃OH and HCONH₂ (as illustrative examples) on the crystalline ice (Section 3.1.1) and on the ASW (Section 3.1.2). The BE values (in kelvin) are reported in each plot.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison between the DFT//HF-3c BEs (in kelvin) computed on the crystalline ice (filled blue circles) and ASW (open circles), respectively, for 20 species studied here: HCl is missing as it dissociates on the ASW (see text).</figcaption></figure>

ASW, as well as the BE associated with each geometry. For methanol, the BE is 8648 K in the crystalline ice, whereas it is 4414 and 10,091 K in the two shown ASW sites. Similarly, for formamide, BE is 9285 K on the crystalline ice and 6639 and 8515 K on the ASW. These two examples show that BEs on ASW can differ more than twice depending on the site and that the value on the crystalline ice can also be substantially different from that on the ASW.

Figure 9 shows the computed BEs, on crystalline ice and ASW, for the studied species. The list of all computed BE values on ASW is reported in Table 2, while Table 3 reports the computed minimum and maximum BE values on ASW and the BEs on the crystalline ice for all the studied species. As already mentioned when presenting the methanol and formamide examples, the amorphous nature of the ice can yield large differences in the calculated BEs with respect to the crystalline values. Figure 9 shows that while the BEs for crystalline versus amorphous ices are very close to each other for H₂, O₂, N₂,
<table><caption>**Table 2**
BE Values (K) Calculated with the DFT//HF-3c Method for Every Case on the Amorphous Slab Model, Where the ZPE Correction Has Not Been Added</caption><thead><tr><th rowspan="2">Species</th><th colspan="8">Amorphous Ice BE Values</th></tr><tr><th>Case 1</th><th>Case 2</th><th>Case 3</th><th>Case 4</th><th>Case 5</th><th>Case 6</th><th>Case 7</th><th>Case 8</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>469</td><td>505</td><td>277</td><td>361</td><td>265</td><td></td><td></td><td></td></tr><tr><td>O<sub>2</sub></td><td>818</td><td>854</td><td>854</td><td>529</td><td>337</td><td></td><td></td><td></td></tr><tr><td>N<sub>2</sub></td><td>1347</td><td>1708</td><td>1311</td><td>1191</td><td>890</td><td></td><td></td><td></td></tr><tr><td>CH<sub>4</sub></td><td>1323</td><td>1960</td><td>1636</td><td>1467</td><td>1070</td><td></td><td></td><td></td></tr><tr><td>CO</td><td>1816</td><td>2189</td><td>1540</td><td>1527</td><td>1299</td><td></td><td></td><td></td></tr><tr><td>CO<sub>2</sub></td><td>2863</td><td>3452</td><td>2538</td><td>2550</td><td>1744</td><td></td><td></td><td></td></tr><tr><td>OCS</td><td>3404</td><td>2971</td><td>2670</td><td>2562</td><td>1780</td><td></td><td></td><td></td></tr><tr><td>HCN</td><td>2923</td><td>5124</td><td>3620</td><td>5136</td><td>3271</td><td>5857</td><td>6146</td><td>7421</td></tr><tr><td>H<sub>2</sub>O</td><td>7156</td><td>5845</td><td>6014</td><td>5689</td><td>4222</td><td></td><td></td><td></td></tr><tr><td>H<sub>2</sub>S</td><td>2814</td><td>3909</td><td>3151</td><td>3560</td><td>2682</td><td></td><td></td><td></td></tr><tr><td>NH<sub>3</sub></td><td>8840</td><td>5268</td><td>6820</td><td>5930</td><td>6579</td><td>6098</td><td>5052</td><td></td></tr><tr><td>CH<sub>3</sub>CN</td><td>8960</td><td>5857</td><td>5617</td><td>5557</td><td>6483</td><td></td><td></td><td></td></tr><tr><td>CH<sub>3</sub>OH</td><td>6531</td><td>4414</td><td>6519</td><td>6362</td><td>5208</td><td>6663</td><td>5509</td><td>10091</td></tr><tr><td>H<sub>2</sub>CO</td><td>3596</td><td>4258</td><td>4174</td><td>4775</td><td>5268</td><td>4594</td><td>3873</td><td>7253</td></tr><tr><td>HCONH<sub>2</sub></td><td>12833</td><td>7481</td><td>10344</td><td>6820</td><td>8467</td><td>7072</td><td>6783</td><td>7313</td></tr><tr><td>HCOOH</td><td>7577</td><td>7409</td><td>8515</td><td>12364</td><td>6302</td><td>7204</td><td>6639</td><td>11354</td></tr><tr><td>OH•</td><td>6230</td><td>1816</td><td>4955</td><td>2754</td><td>5076</td><td></td><td></td><td></td></tr><tr><td>HCO•</td><td>2694</td><td>2057</td><td>1540</td><td>3608</td><td>1672</td><td></td><td></td><td></td></tr><tr><td>CH<sub>3</sub>•</td><td>1708</td><td>1936</td><td>1299</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>NH<sub>2</sub>•</td><td>4354</td><td>3716</td><td>4402</td><td>3368</td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>**Table 3**
Summary of Our Computed BEs and Comparison with Data from the Literature</caption><thead><tr><th rowspan="3">Species</th><th colspan="4">BEs from This Work</th><th colspan="4">BEs from Literature</th></tr><tr><th rowspan="2">Crystalline Ice<br>BE(0) disp</th><th colspan="2">ASW</th><th colspan="2">Computed</th><th colspan="2">Databases</th><th colspan="2">Experiments</th></tr><tr><th>Min</th><th>Max</th><th>Das<sup>(a)</sup></th><th>Wakelam<sup>(b)</sup></th><th>UMIST<sup>(c)</sup></th><th>KIDA<sup>(d)</sup></th><th>Penteado<sup>(e)</sup></th><th>Others</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>790</td><td>226</td><td>431</td><td>545</td><td>800</td><td>430</td><td>440</td><td>480 ± 10</td><td>322–505<sup>(f)</sup></td></tr><tr><td>O<sub>2</sub></td><td>677</td><td>287</td><td>729</td><td>1352</td><td>1000</td><td>1000</td><td>1200</td><td>914–1161</td><td>920–1520<sup>(f)</sup></td></tr><tr><td>N<sub>2</sub></td><td>1242</td><td>760</td><td>1458</td><td>1161</td><td>1100</td><td>790</td><td>1100</td><td>1200</td><td>790–1320<sup>(f1)</sup></td></tr><tr><td>CH<sub>4</sub></td><td>1633</td><td>914</td><td>1674</td><td>2321</td><td>800</td><td>1090</td><td>960</td><td>1370</td><td>900–1800<sup>(f1)</sup></td></tr><tr><td>CO</td><td>1663</td><td>1109</td><td>1869</td><td>1292</td><td>1300</td><td>1150</td><td>1300</td><td>863–1420</td><td>960–1947<sup>(f1–h)</sup></td></tr><tr><td>CO<sub>2</sub></td><td>2568</td><td>1489</td><td>2948</td><td>2352</td><td>3100</td><td>2990</td><td>2600</td><td>2236–2346</td><td></td></tr><tr><td>OCS</td><td>2722</td><td>1520</td><td>2907</td><td>1808</td><td>2100</td><td>2888</td><td>2400</td><td>2325<sup>a</sup></td><td>2430<sup>(i)</sup></td></tr><tr><td>HCl</td><td>5557</td><td>(l)</td><td>(l)</td><td>4104</td><td>4800</td><td></td><td>5172</td><td></td><td>5172<sup>(m)</sup></td></tr><tr><td>HCN</td><td>6392</td><td>2496</td><td>6337</td><td>2352</td><td>3500</td><td>2050</td><td>3700</td><td></td><td></td></tr><tr><td>H<sub>2</sub>O</td><td>7200</td><td>3605</td><td>6111</td><td>4166</td><td>4600</td><td>4800</td><td>5600</td><td>4815–5930</td><td></td></tr><tr><td>H<sub>2</sub>S</td><td>4468</td><td>2291</td><td>3338</td><td>3232</td><td>2500–2900</td><td>2743</td><td>2700</td><td>2296<sup>a</sup></td><td></td></tr><tr><td>NH<sub>3</sub></td><td>6235</td><td>4314</td><td>7549</td><td>5163</td><td>5600</td><td>5534</td><td>5500</td><td>2715<sup>a</sup></td><td></td></tr><tr><td>CH<sub>3</sub>CN</td><td>5906</td><td>4745</td><td>7652</td><td>3786</td><td>4300</td><td>4680</td><td>4680</td><td>3790<sup>a</sup></td><td></td></tr><tr><td>CH<sub>3</sub>OH</td><td>7385</td><td>3770</td><td>8618</td><td>4511</td><td>4500–5100</td><td>4930</td><td>5000</td><td>3820<sup>a</sup></td><td>3700–5410<sup>(n,o,p)</sup></td></tr><tr><td>H<sub>2</sub>CO</td><td>5187</td><td>3071</td><td>6194</td><td>3242</td><td>5100</td><td>2050</td><td>4500</td><td>3260 ± 60</td><td></td></tr><tr><td>HCONH<sub>2</sub></td><td>8104</td><td>5793</td><td>10960</td><td></td><td>6300</td><td>5556</td><td></td><td></td><td>7460–9380<sup>(q)</sup></td></tr><tr><td>HCOOH</td><td>7991</td><td>5382</td><td>10559</td><td>3483</td><td>5000</td><td>5570</td><td>4532<sup>a</sup></td><td></td><td></td></tr><tr><td>OH•</td><td>5588</td><td>1551</td><td>5321</td><td>3183</td><td>3300–5300</td><td>2850</td><td>4600</td><td>1656–4760</td><td></td></tr><tr><td>HCO•</td><td>2968</td><td>1315</td><td>3081</td><td>1857</td><td>2300–2700</td><td>1600</td><td>2400</td><td></td><td></td></tr><tr><td>CH<sub>3</sub>•</td><td>2188</td><td>1109</td><td>1654</td><td>1322</td><td>2500</td><td>1175</td><td>1600</td><td></td><td></td></tr><tr><td>NH<sub>2</sub>•</td><td>5156</td><td>2876</td><td>4459</td><td>3240</td><td>2800–4500</td><td>3956</td><td>3200</td><td></td><td></td></tr></tbody></table>

**Notes.** Column (1) reports the species, Columns (2)–(4) the BEs computed in the present work and corrected for the ZPE, Columns (5) and (6) the values obtained via calculations from other authors, Columns (7) and (8) the values in the two astrochemical databases KIDA and UMIST (see text), and Columns (9) and (10) the values measured in different experiments. Units are in K, and the references are listed in the notes below. References: (a) Das et al. 2018; (b) Wakelam et al. (2017); (c) McElroy et al. (2013); (d) Wakelam et al. (2015); (e)Penteado et al. (2017); (f) He et al. (2016), note that (f1) refers to porous ice; (g) Raut et al. (2007); (h) Smith et al. (2016); (i) Ward et al. (2012); (l) HCl molecule dissociate; (m) Olanrewaju et al. (2011); (n) Minissale et al. (2016); (o) Martín-Doménech et al. (2014); (p) Bahr et al. (2008); (q) Chaabouni et al. (2018), note that the BE refers to the silicate substrate because it is larger than that of water ice.

<footnote>a Results estimated from the work of Collings et al. (2004), reported in Table 2 of Penteado et al. (2017).</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Spider graphs of the DFT/HF-3c BE values (in kelvin) calculated on the same eight adsorption sites of the ASW for the H₂CO (red), CH₃OH (black), HCONH₂ (green), and HCOOH (blue) molecules. The BE value scale goes from 0 K (center of the graph) to 14,433 K (vertices of the polygon) in steps of 4811 K. Labeling of dH and dO sites is referring to Figure 3.</figcaption></figure>

CH₄, CO, CO₂, and OCS, the ones computed for the remaining molecules for the crystalline ice fall in the highest range of the distribution of the amorphous BE values. This behavior can be explained considering the smaller distortion energy cost upon adsorption for the crystalline ice compared to the amorphous one. The different local environment provided by crystalline versus amorphous ices is also the reason for HCl being molecularly adsorbed at the crystalline ice while becoming dissociated at the amorphous one. Further details about the case of HCl are reported in Appendix A.2. This probably will not occur for HF, not considered here, which is expected to be molecularly adsorbed on both ices owing to its higher bond strength compared to HCl. Nevertheless, as we did not explore exhaustively all possible configurations of the adsorbates at the amorphous surfaces, we cannot exclude that some even more/less energetic binding cases remain to be discovered.

Some adsorbates show similar trends in the BEs, despite their different chemical nature. This is shown in Figure 10, in which we plot the BE values on ASW for four molecules that have been adsorbed at the same adsorption sites: formaldehyde, formic acid, formamide, and methanol. The BE distributions for the H₂CO and HCOOH are very similar (in their relative values), and those for CH₃OH and HCONH₂ show some similarities, despite the large difference in the chemical functionality.

# 4. Discussion

A first rather expected result of our computations is that the BE of a species on the ASW is not a single value: depending on the species and the site where it lands, the BE can largely differ, even by more than a factor of two (Table 3). This has already been discussed in the literature, for instance, for H adsorption on both crystalline and amorphous ice models (Ásgeirsson et al. 2017). This has important consequences both when comparing the newly computed BE(0)s with those in the literature Section 4.1 and for the astrophysical implications in

Section 4.2. We will discuss these two aspects separately in the next two sections.

## 4.1. Comparison BE Values in the Literature

Being such a critical parameter, BEs have been studied from an experimental and theoretical point of view. In this section, we will compare our newly computed values with those in the literature, separating the discussion for the experimental and theoretical values, respectively. We will then also comment on the values available in the databases that are used in many astrochemical models.

### 4.1.1. Comparison with Experimental Values

In the present computer simulation we have computed the BE released when a species is adsorbed on the surfaces of the ice models (either crystalline or amorphous) at very low adsorbate coverage $\theta$ ($\theta \to 0$). The correct comparison with experiments would therefore be with microcalorimetric measurements at the zero-limit adsorbate coverage. In astrochemical laboratories, TPD is, instead, the method of choice and is related to the desorption activation energy (DAE). DAE derives indirectly from the TPD peaks through Readhead's method (Redhead 1962), or more sophisticated techniques. TPD usually starts from an ice surface hosting a whole monolayer of the adsorbate and therefore depends also on $\theta$ (He et al. 2016), rendering the comparison with the theoretical BE not straightforward (King 1975). Ice restructuring processes may also affect the final DAE. Sometimes, TPD experiments only provide desorption temperature peaks $T_{des}$, without working out the DAE. This is the case of the fundamental work by Collings et al. (2004). Therefore, BEs reported in the review by Penteado et al. (2017) relative to the Collings et al. data (see Table 3) were computed through the approximate formula: $BE(X) = [T_{des}(X)/T_{des}(H_2O)] BE(H_2O)$, in which $T_{des}(X)$ is the desorption temperature of the X species contrasted with that of water $T_{des}(H_2O)$ to arrive at the corresponding $BE(X)$ by assuming that of water to be 4800 K. For the above reasons, a one-by-one comparison between experiment and modeling is outside the scope of the present paper.

Following the above warnings, we can now analyze Table 3 reporting the recent compilation by Penteado et al. (2017) (vide supra) plus the values that appeared in the literature after that compilation. We start with the cases of two measurements carried out by He et al. (2016) on porous amorphous ice surfaces, for N₂ and CO. Table 3 shows two values reported by He et al. for the two extreme cases of when the ice is completely covered by the species (the smaller value) and when, on the contrary, it is less than a monolayer (the largest value), which is the one to compare with our computed values. Our BEs on amorphous ice models are in reasonable agreement with those measured by He et al. for CO (1109–1869 K vs. 1940 K) and on the lower side for N₂ (760–1455 K vs. 1800 K). It is worth noting that the comparison is much better when referring to the nonporous amorphous ice measurements by the same authors: for CO and N₂ the measured BE values are 1600 and 1320 K, indeed well bracketed by our BE(0). Data from Penteado et al. (2017) extracted from the TPD of Collings et al. (2004) for NH₃ adsorbed on the ice layer gave a BE of 2715 K. This value is, however, identical to that from the TPD of NH₃ adsorbed on the gold surface (no water ice), proving that that BE is relative to the NH₃/NH₃ lateral interaction within the
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison of the computed zero-point energy corrected BE(0)s for the amorphous ice model with respect to those by Das et al., by Wakelam et al., and reported in the KIDA and UMIST databases (McElroy et al. 2013; Wakelam et al. 2015, 2017; Das et al. 2018).</figcaption></figure>

adsorbed NH₃ multilayer and not due to the interaction with the ice surface. For reasons explained in Collings et al. (2004), the BE of NH₃ on ice is assumed of the same order as that of water, i.e., around 4800 K, in better agreement with the UMIST value of 5500 K (McElroy et al. 2013). Indeed, our data (see Table 3, ZPE corrected) of 4300–7500 K bracket the experimental ones. The computed highest values emphasize the H-bond acceptor capability of NH₃ occurring on a few specific sites characterized by very high electrostatic potential, only important for very low NH₃ coverage, not easily accessible in the TPD experiments. For the H₂O case, the computed BE(0)s (Table 3) for the amorphous ice are in the 3605–6111 K range, reasonably bracketing the experimental one of 4815–5930 K.

In general, the comparison of our BE values computed on the ASW with those measured by the various experiments reported in Table 3 shows an excellent agreement, when considering the ranges in our values and the ranges in the values of the experiments. Only one species seems to have relatively different computed and measured BEs: O₂. For O₂, experiments tend to provide larger values with respect to what we computed (our largest value is 729 K, while the lowest measured value is 914 K). For many other species, except H₂, our computed lowest BE(0)s are within the range of the measured ones, but we predict sites

where BE(0)s are larger, which may have important astrophysical implications (Section 4.2). Finally, for H₂ we predict sites where the BE(0) is (slightly) lower than the measured ones.

### 4.1.2. *Comparison with Computed Values*

In the literature, there are two works that reported computations of BEs for a large set of molecules, those by Wakelam et al. (2017) and those by Das et al. (2018). The former carried out computations considering only one water molecule, whereas the latter considered a cluster of up to six water molecules.

The first aspect to notice is, therefore, that neither of these two studies can, by definition, reproduce the strong adsorption sites that we have in our ASW model. Indeed, only the adoption of more realistic and periodically extended ice models allows us to fully consider the hydrogen bond cooperativity, which will enhance the strength of the interaction with adsorbates at the terminal dH atoms exposed at the surface. This important effect is entirely missed by the two above-mentioned works. It is not surprising, then, that our crystalline and ASW BEs differ, sometimes substantially, from the Wakelam et al. and Das et al. values (as, by the way, they differ between themselves as well). This is clearly shown in Figure 11, where we report the comparison of our computed values with those by Wakelam et al. and Das et al., respectively. In general,
both work values tend to lay in the low end of our computed BEs. As extreme examples, our ASW BEs are larger for CH₃CN and HCOOH. The inverse effect is observed for the smallest studied species: our BEs are smaller than those computed by Wakelam et al. and Das et al. for H₂ and O₂.

### 4.1.3. Comparison with Values in Astrochemical Databases

Two databases list the BEs of the species used by the astrochemical models: the Kinetic Database for Astrochemistry (KIDA, http://kida.astrophy.u-bordeaux.fr/; Wakelam et al. 2015) and UMIST (http://udfa.ajmarkwick.net/index.php; McElroy et al. 2013). The comparison between our newly computed values and those reported in the two databases is shown in Figure 11. The general remarks that we wrote for the comparison with the literature experimental and theoretical values (Sections 4.1.1 and 4.1.2) roughly apply here: the databases quote BE values in the low end of ours. This is not surprising, as the databases are compiled based on the experimental and theoretical values in the literature. We just want to emphasize here, once again, that the sites with large BEs are lacking, and this may have important consequences in the astrochemical model predictions.

## 4.2. Astrophysical Implications

BEs enter in two hugely important ways in the chemical composition of interstellar objects/clouds: (i) they determine at what dust temperature the frozen species sublimate, and (ii) they determine at what rate the species can diffuse in the ice, as the diffusion energy is a fraction of the BE species. Both processes are mathematically expressed by an exponential containing the BE. Therefore, even relatively small variations of the BE can cause huge differences in the species abundances in the gas phase and on the grain surfaces, where they can react with other species.

In this context, probably the most important astrophysical implication of the present study is that in our ASW model (which is likely the best description to represent the interstellar amorphous ice so far available in the context of the BE estimates) a species does not have a single value, but a range of values that depend on the species itself and the site where it is bound. The range can spread by more than a factor of two: this obviously can have a huge impact on the modeling and, consequently, our understanding of the interstellar chemical evolution.

### 4.2.1. Impact of Multiple BEs in Astrochemical Modeling

To give a practical example of the impact on the gaseous abundance, we built a toy model for the interstellar ice and simulated the desorption rate of the ice as a function of the temperature. Our scope here is not to compare the toy model predictions with astronomical observations or laboratory experiments: we only mean to show how multiple BEs (we used the electronic BEs) would lead, in principle, to a different behavior of the ice sublimation process. Therefore, we developed a toy model that does not contain diffusion or reaction processes on the ice surface or rearrangement of the ice during the ice heating, but only a layered structure with two species, specifically water and methanol, where molecules have the range of BEs calculated in Section 3. We then show how the multiple BEs affect the temperature at which peaks of desorption appear, considering that only species at the surface

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Desorption rate of methanol (red curves) and water (blue curves) as a function of the temperature. In these computations, the ice is assumed to be composed of 10 layers: the bottom five layers contain water, while the top five contain 80% water and 20% methanol. The ice is assumed to be at 10 K at the beginning of the simulation, and it reaches 400 K in $10^5$ yr. The curves refer to the case when a single BE (from the KIDA database) is considered for water and methanol, respectively (dashed line), and when the multiple BEs of this work are considered (solid line). The bottom panel shows the methanol and water desorption rates normalized to 1.</figcaption></figure>

of the ice, namely, exposed to the void, can sublimate (and not the entire bulk). In this toy model, we considered 10 layers for an icy grain mantle, where the bottom five are made entirely of water and each of the top five layers contains 80% water and 20% methanol. The methanol molecules with different BEs are distributed randomly on each layer, with the same proportion of BE sites. Looking at the BE values computed with the DFT//HF-3c method on the ASW model, methanol has eight BE values (4414, 5208, 5509, 6362, 6519, 6531, 6663, and 10,091 K), and at each layer there will be 12.5% methanol molecules with each of the eight BEs. The same applies for water molecules, for which we have computed five possible BEs (4222, 5689, 5845, 6014, and 7156 K), with 20% water molecules with each BE. In this model, only molecules of the layers in contact with the void can evaporate: for example, methanol molecules can be trapped if they have water molecules with larger BEs on top of them.

We start with an ice temperature of 10 K, and at the end it reaches 400 K in $10^5$ yr, to simulate the heating of a collapsing solar-like protostar. The plot of the desorption rates is shown in Figure 12, where we also show them assuming the BE values from the KIDA database for methanol and water, respectively. First, when only the KIDA values for BE are assumed, water molecules desorb at about 110 K; methanol has two peaks of desorption rates, the first at about 95 K, corresponding to the desorption of the methanol molecules not trapped by the water molecules, and the second peak at about 110 K, when all water molecules desorb so that no methanol molecules are trapped. Note that the desorption of the water molecules of the bottom layers arrives at a slightly larger temperature.
Not surprisingly, the introduction of multiple BEs produces multiple peaks of desorption, for both water and methanol. Figure 12 shows that the water desorption rate has a small peak at ~75 K, a larger one at ~120 K, then another at ~140 K, and, finally, a last peak at ~190 K. Methanol starts to desorb at 80 K, the bulk is desorbed between 120 and 140 K, and a last peak is seen around 190 K. We emphasize again that this is a toy model meant to show the potential impact of the new BEs on the astrochemical modeling. The details will depend on the real structure of the water ice and how molecules are distributed on the icy mantles. They will determine how many sites have a certain BE value and how molecules are trapped in the ice. As a very general remark, we can conclude that species can be in the gas phase at lower and higher dust temperatures than if one only considers a single BE.

### 4.2.2. Looking Forward: Implementation of Multiple BEs in Astrochemical Models

Our toy model introduced in the previous subsection shows the importance of considering multiple BEs for each species in the astrochemical models to have more realistic predictions. In this work, we provide the possible BEs for 21 species (Table 3). Very likely, they cover most of the possibilities, as they span a large range of H-bonds within the water molecules of the ASW. However, from a computational standpoint, such an adsorption variability has to be fully explored, in which plotting the different calculated BEs in histograms is useful to provide insights on the shape of the BE distribution (Song & Kästner 2017). Moreover, in order to build a reliable astrochemical model, one would also need to know the relative probability for each BE, and our present study is unable to provide sensible numbers. For that, a statistical study on an ASW model that is much larger than the one used here is necessary. This is a step that we indeed plan to take in future studies. Meanwhile, we adopted a distribution in which we assign an equal fraction of molecules to each BE. If one looks carefully at the distribution of the BEs for each molecule, they are not uniformly distributed but peak around some values: for example, methanol has a peak around 6000 K and an extreme value of 10,091 K only, so that, very likely, this site will be less populated than the sites around 6000 K, as shown by Figure 10. Yet, considering even a smaller fraction of these extreme values may have important consequences, for example, in the so-called snow lines of protoplanetary disks, or even on the observed abundances toward hot cores and hot corinos, or, finally, toward prestellar cores.

### 4.2.3. Comments on N₂, CO, and HCl

Finally, we would like to comment on three species of the studied list, N₂, CO, and hydrogen chloride (HCl).

**N₂ and CO:** Our computations show that the BE of CO is definitively larger than that of N₂, against the values that are in the astrochemical databases (see Table 3): on average, our computed BEs differ by about 400 K, whereas in the databases the difference is 200 and 360 K in KIDA and UMIST, respectively. This difference very likely can explain why observations detect N₂H⁺, which is formed in the gas phase from N₂, where CO is already frozen on the grain mantles (e.g., Bergin et al. 2002; Tafalla et al. 2004; Redaelli et al. 2019), a debate that has been going on for almost two decades (e.g., Öberg et al. 2009; Pagani et al. 2012). In order to quantify the

effect, a specifically focused modeling will be necessary, which is beyond the scope of the present work. Here we want simply to alert that the new BEs might explain some long-standing mysteries. Another comment regards the difference in the BEs on crystalline surfaces and ASW. Again, the CO BE is about 500 K larger than that of N₂, and both are larger by about 300 K than those on the ASW, a difference that also has an impact on the snow lines of these two species in protoplanetary disks, where crystalline water ices have been detected (Terada & Tokunaga 2012).

**HCl:** Astrochemical models predict that HCl is the reservoir of Cl in molecular gas (e.g., Schilke et al. 1995; Neufeld & Wolfire 2009; Acharyya & Herbst 2017). However, all the observations carried out so far have found that only a tiny fraction of Cl is in the gaseous HCl, even in sources where all the grain mantles are supposed to be completely sublimated (Peng et al. 2010; Codella et al. 2011; Kama et al. 2015). One possible explanation is that HCl, once formed in the gas phase, is adsorbed on the grain icy mantles and dissociates, as shown by our calculations on the ASW model and also by previous calculations on the crystalline P-ice model (Casassa 2000) and for more sophisticated proton-disordered crystalline ice models (Svanberg et al. 2000). It is a matter to be studied whether the sublimation of the water, when the dust reaches about 100–120 K, would also provide a reactive channel transforming the Cl anion in the neutral atom, the latter obviously unobservable. This would help in solving the mystery of HCl not being observed in gas phase. Furthermore, if that were the case, the population of the chemically reactive atomic Cl would be increased, with an important role in the gas-phase chemistry (see, e.g., Balucani et al. 2015; Skouteris et al. 2018).

## 5. Conclusions

In this work, we present both a new computational approach and realistic models for crystalline and amorphous water ice to be used to address an important topic in astrochemistry: the BEs of molecules on interstellar ice surfaces. We simulated such surfaces by means of two (antipodal) models, in both cases adopting a periodic approach: a crystalline and an amorphous 2D slab model. We relied on DFT calculations, using the B3LYP-D3 and M06-2X widely used functionals. This approach was further validated by an ONIOM-like correction at the CCSD(T) level. Results from this combined procedure confirm the validity of the BEs computed with the adopted DFT functionals. The reliability of a cost-effective HF-3c method adopted to optimize the structures at the amorphous ice surface sites was proved by comparing the BEs computed at the crystalline ice surface at the DFT//DFT and DFT//HF-3c levels, which were found to be in very good agreement.

On both ice surface models, we simulated the structure and adsorption energetic features of 21 molecules, including 4 radicals, representative of the most abundant species of the dense ISM. A main conclusion is that the crystalline surfaces only show very limited variability in the adsorption sites, whereas the amorphous surfaces provide a wide variety of adsorption binding sites, resulting in a distribution of the computed BE. Furthermore, BE values at crystalline ice surface are in general higher than those computed at the amorphous ice surfaces. This is largely due to the smaller geometry relaxation cost upon adsorption compared to the amorphous cases, imposed by the tighter network of interactions of the denser crystalline ice over the amorphous ice.
Finally, the BEs obtained by the present computations were compared with literature data, from both experimental and computational works, as well as those on the public astrochemical databases KIDA and UMIST. In general, our BEs agree relatively well with those measured in the laboratory, with the exception of O₂ and, to a lesser extent, H₂. On the contrary, previous computations of BEs, which considered a very small number of (≤6) water molecules, provide generally lower values with respect to our new computations and, with no surprise, miss the fact that BEs have a spread of values that depend on the position of the molecule on the ice. Since the two astrochemical databases mentioned above are based on the literature data, our BEs differ, sometimes substantially, from those quoted and do not report multiple BE values.

We discussed some astrophysical implications, showing that the multiple computed BEs give rise to a complex process of interstellar iced mantle desorption, with multiple peaks as a function of the temperature that depends (also) on the ice structure. Our new computations do not allow us to estimate how the BEs are distributed for each molecule, as only a statistical study is necessary for that. The new (multiple) BEs of N₂ and CO might explain why N₂H⁺ depletes later than CO in prestellar cores, while the relatively low abundance of HCl, observed in protostellar sources, could be due to the fact that it dissociates into the water ice, as shown by our calculations.

Finally, the present study shows the importance of theoretical calculations of BEs on as realistic as possible ice surfaces. This first study of 21 molecules needs to be extended to the hundreds of molecules that are included in the astrochemical models to have a better understanding of the astrochemical evolution of the ISM.

A part of the computational results were from the SF Master thesis “Ab initio quantum mechanical study of the interaction of astrochemical relevant molecules with interstellar ice models,” Dipartimento di Chimica, University of Torino, Torino, 2018. S.F., L.Z., and P.U. acknowledge financial support from the Italian MIUR (Ministero dell'Istruzione, dell'Università e della Ricerca) and from Scuola Normale Superiore (project PRIN 2015, STARS in the CAOS—Simulation Tools for Astrochemical Reactivity and Spectroscopy in the Cyberinfrastructure for Astrochemical Organic Species, cod. 2015F59J3R). The Italian CINECA consortium is also acknowledged for the provision of supercomputing time for part of this project. A.R. is indebted to the “Ramón y Cajal” program. MINECO (project CTQ2017-89132-P) and DIUE (project 2017SGR1323) are acknowledged. This project has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 811312 for the project “Astro-Chemical Origins” (ACO) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program, for the Project “The Dawn of Organic Chemistry” (DOC), grant agreement No. 741002. Finally, we wish to acknowledge the extremely useful discussions with Prof. Gretobape.

## Appendix A
Computational Details

In CRYSTAL17, the multielectron wave function is built as a Slater determinant of crystalline/molecular orbitals, which are linear combinations of localized functions on the different atoms of the structure that are called atomic orbitals (AOs). In a similar manner, the AOs are constructed by linear combinations of localized Gaussian functions that form a basis set. The basis set employed for this work is an Ahlrichts-TVZ (Schäfer et al. 1992), added with polarization functions.

### A.1. BEs, Counterpoise, and Zero-point Energy Corrections

In a periodic treatment of surface adsorption phenomena one of the most relevant energy values, useful to describe the interacting system, is the BE, which is related to the interaction energy $ΔE$, so that

$$BE = -\Delta E. \quad (A1)$$

The BE per unit cell per adsorbate molecule BE is a positive quantity (for a bounded adsorbate) defined as

$$BE = [E_m(M // M) + E(S // S)] - E(SM // SM), \quad (A2)$$

where $E(SM // SM)$ is the energy of a fully relaxed unitary cell containing the surface slab S in interaction with the adsorbate molecules M, $E(S // S)$ is the energy of a fully relaxed unitary cell containing the slab alone, and $E_m(M // M)$ is the molecular energy of the free fully optimized adsorbate molecule (the symbol following the double slash identifies the geometry at which the energy, E, is calculated)

$$BE = BE^* - \delta E_S - \delta E_M \quad (A3)$$

$$\delta E_S = E(S // SM) - E(S // S) \quad (A4)$$

$$\delta E_M = E(M // SM) - E_m(M // M) \quad (A5)$$

$$BE^* = [E(S // SM) + E(M // SM)] - E(SM // SM), \quad (A6)$$

in which $\delta E_S$ is the deformation energy of the surface ($\delta E_S > 0$), whereas $\delta E_M$ ($= \Delta E_M + \Delta E_L$) counts both the deformation energy of the molecule ($\Delta E_M$) and the lateral intermolecular interactions ($\Delta E_L$) between the infinite molecule images in the same spatial configuration occurring in the SM periodic system. The pure molecule’s deformation energy can be computed as

$$\Delta E_M = E_m(M // SM) - E_m(M // M), \quad (A7)$$

in which $E_m(M // SM)$ is the molecular energy of the molecule frozen at the geometry occurring on the surface and $E_m(M // M)$ is the molecular energy of a fully optimized free molecule, so that $\Delta E_M > 0$. The lateral intermolecular interactions, $\Delta E_L$, are defined as

$$\Delta E_L = E(M // SM) - E_m(M // SM) \quad (A8)$$

and can be either positive (repulsion) or negative (attraction). With those positions, the $BE^*$ interaction energy is then deformation and lateral interaction free, being the result of energy differences between periodic calculations carried out at the geometry of the SM system. For instance, $E(M // SM)$ is the energy of the unit cell of a crystal containing only a molecule in
<figure><img src="image_14.png" /><figcaption>Figure 13. Linear regression between BE and BE(0) (zero-point-corrected) values computed for the considered adsorbates on the crystalline ice.</figcaption></figure>

the same geometry assumed in the SM system. The above BE definition can be easily recast to include the BSSE correction, using the same counterpoise method adopted for intermolecular complexes (Boys & Bernardi 1970; Davidson & Feller 1986). The definition of the BSSE-corrected interaction energy BEC is then

$$BE^C = BE^{*C} - \delta E_S - \Delta E_M - \Delta E_L^C \quad (A9)$$

$$BE^{*C} = [E(S[M]//SM) + E([S]M//SM)] - E(SM//SM) \quad (A10)$$

$$BSSE = BE - BE^C, \quad (A11)$$

in which E(S[M]//SM) and E([S]M//SM) are the energy of the slab plus the ghost functions of the molecules and the energy of the infinite replica of molecules with the ghost functions of the underneath slab, respectively. Because the variational theorem ensures that BEC < BE*, it immediately results that BSSE > 0.

Each of the terms of Equation (A2) can be corrected by the zero-point vibrational contribution (in the harmonic approximation), ZPE, therefore arriving at the definition of the zero-point correct BE, BE(0), as

$$BE(0) = [E_m(M // M) + E(S // S)] - E(SM // SM) + [ZPE(M) + ZPE(S) - ZPE(SM)] \quad (A12)$$

$$BE(0) = BE - \Delta ZPE \quad (A13)$$

$$\Delta ZPE = ZPE(SM) - ZPE(M) - ZPE(S). \quad (A14)$$

In this work the ΔZPE has been evaluated for the crystalline ice cases only, and then the scaling factor found in the linear fit of Figure 13 was adopted to correct the BEs for the amorphous ice model for the reason described in Section 2.1.

<figure><img src="image_15.png" /><figcaption>Figure 14. Enlarged views of HCl adsorbed on a crystalline and amorphous ice models.</figcaption></figure>

## A.2. Details on HCl Adsorbed on Crystalline and Amorphous Ices

Our calculations showed that HCl remains molecularly adsorbed at the crystalline ice surface. At the amorphous one it spontaneously deprotonates, making a locally stable ion pair (Cl-/H3O+) (see the figure below). The reason is that for the surface selected to represent the crystalline ice, the ion pair cannot be stabilized by a large enough network of H-bonds owing to the rigidity of the structure. On the contrary, at amorphous surface, Cl- is engaged in three H-bonds, while the hydronium ion remains embedded in the H-bonds provided by the ice surface. A good solvation of the ion pair is the key determining the final fate of HCl at the ice surfaces (molecular or ion pair) as already pointed out many years ago by Novoa and Sosa for a small water cluster hosting HCl (Novoa & Sosa 1995). While deprotonation of HCl cannot be excluded for specific crystalline ice surface sites as simulated by Svanberg et al. (2000) when enough dH can stabilize the anion, deprotonation would be much more common at amorphous ice surface or locally distorted crystalline ones, due to the presence of favorable local environment (see Figure 14). Dissociative adsorption was found early by Horn et al. (1992) through RAIRS spectroscopy of DCl adsorption on thin D2O film, and by Olanrewaju et al. (2011) through thermal and electron-stimulated desorption experiments.

## A.3. CRYSTAL17 Computational Parameters

In order to optimize the values of the shrinking factor, the tolerances on integrals of the integration grid (SHRINK, TOLINTEG, and GRID parameters in the code; Dovesi et al. 2018) as described in Section 2 of the article, NH3 has been adopted as a test case. Geometry optimizations for the adsorption process on the crystalline slab model have been run with the previously described computational scheme, varying one by one the aforementioned parameters, with convergence on the pure ΔE as defined by Equation (A1) being the threshold. For this benchmark calculation, BSSE, distortion, and lateral interaction contributions have not been computed. Results are resumed in what follows.

From Table 4 it is clear that ΔE is practically unaffected by the SHRINK value; thus, SHRINK = (2 2) has been used for all calculations, allowing the saving of computational time.
<table><caption>**Table 4**
Optimization of the SHRINK Parameter</caption><thead><tr><th>SHRINK</th><th>k Points in BZ</th><th>ΔE (K)</th></tr></thead><tbody><tr><td>2 2</td><td>4</td><td>7873</td></tr><tr><td>3 3</td><td>5</td><td>7889</td></tr><tr><td>4 4</td><td>10</td><td>7889</td></tr><tr><td>5 5</td><td>13</td><td>7890</td></tr><tr><td>6 6</td><td>20</td><td>7890</td></tr><tr><td>7 7</td><td>25</td><td>7890</td></tr><tr><td>8 8</td><td>34</td><td>7890</td></tr></tbody></table>

**Note.** Tolerances of integrals (TOLINTEG) and integration grid (XLGRID) are fixed to default values.

<table><caption>**Table 5**
Optimization of the TOLINTEG Parameter (SHRINK = 2 2 and XLGRID)</caption><thead><tr><th>TOLINTEG</th><th>ΔE (K)</th></tr></thead><tbody><tr><td>6 6 6 6 12</td><td>7889</td></tr><tr><td>7 7 7 7 14</td><td>7906</td></tr><tr><td>8 8 8 8 16</td><td>8063</td></tr><tr><td>9 9 9 9 18</td><td>7981</td></tr><tr><td>10 10 10 10 20</td><td>7890</td></tr></tbody></table>

As for the SHRINK parameter, the variations of TOLINTEG values (see Table 5) do not practically affect the final ΔE value. Consequently, its values have been set equal to 7, 7, 7, 7, and 14.

Again, the ΔE is practically unaffected by changing the adopted grid (see Table 6); therefore, XLGRID was selected for all the remaining calculations. Resuming, every other calculation has been carried out with SHRINK = (2 2), TOLINTEG = (7 7 7 7 14), and XLGRID as computational parameters.

### A.4. Description of Dispersive Forces

We optimized the geometry of two parallel benzene rings that interact with each other just because of dispersive forces

<figure><img src="image_16.png" /></figure>

**Figure 15.** Initial and optimized geometry of two parallel benzene rings obtained with B3LYP (without dispersion corrections), B3LYP-D3 and M06-2X functionals. Distances between the planes of the two rings are reported in angstroms.

<table><caption>**Table 6**
Optimization of the Integration Grid Parameter (SHRINK = 2 2 and TOLINTEG = 7 7 7 7 14)</caption><thead><tr><th>GRID</th><th>ΔE (K)</th></tr></thead><tbody><tr><td>XLGRID</td><td>7951</td></tr><tr><td>XXLGRID</td><td>7951</td></tr></tbody></table>

given the apolar nature of the molecule. The results of Figure 15 clearly show that the B3LYP functional needs some a posteriori correction (like Grimme's D3) for dispersive forces in order to describe correctly the interaction (Grimme et al. 2010), whereas the M06-2X functional is able to describe it correctly without any a posteriori correction (Zhao & Truhlar 2008).
<table><caption>**Table 7**
Summary of the Components Used for the Calculation of the ONIOM2 Correction to the BEs</caption><thead><tr><th>Species</th><th>BE(CCSD(T), μ-nH<sub>2</sub>O) Extrapolated</th><th>BE(DFT, μ-nH<sub>2</sub>O)</th><th>BE(DFT, ice)</th><th>BE(ONIOM2)</th><th>| Diff |</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>843</td><td>796</td><td>1192</td><td>1240</td><td>47</td></tr><tr><td>CO</td><td>1127</td><td>1165</td><td>2356</td><td>2318</td><td>38</td></tr><tr><td>CO<sub>2</sub></td><td>1939</td><td>1817</td><td>3442</td><td>3564</td><td>121</td></tr><tr><td>OCS</td><td>1354</td><td>1138</td><td>3478</td><td>3695</td><td>216</td></tr><tr><td>H<sub>2</sub>S</td><td>2443</td><td>2875</td><td>5679</td><td>5248</td><td>432</td></tr><tr><td>HCN</td><td>3300</td><td>3453</td><td>5795</td><td>5642</td><td>153</td></tr><tr><td>H<sub>2</sub>CO</td><td>3199</td><td>3114</td><td>6490</td><td>6575</td><td>85</td></tr><tr><td>HCl</td><td>3993</td><td>4475</td><td>6501</td><td>6019</td><td>482</td></tr><tr><td>NH<sub>3</sub></td><td>4246</td><td>4389</td><td>7376</td><td>7233</td><td>143</td></tr><tr><td>CH<sub>3</sub>OH</td><td>5196</td><td>5275</td><td>8681</td><td>8602</td><td>79</td></tr><tr><td>HCOOH</td><td>6263</td><td>6223</td><td>9522</td><td>9562</td><td>40</td></tr><tr><td>HCONH<sub>2</sub></td><td>5080</td><td>5056</td><td>9614</td><td>9638</td><td>24</td></tr><tr><td>OH•</td><td>4492</td><td>4502</td><td>6542</td><td>6532</td><td>10</td></tr><tr><td>NH<sub>2</sub>•</td><td>3021</td><td>2788</td><td>6043</td><td>6276</td><td>233</td></tr><tr><td>HCO•</td><td>2339</td><td>2088</td><td>3474</td><td>3725</td><td>251</td></tr></tbody></table>

Note. The last column shows the absolute value of the difference between the BE(DFT, μ-2H₂O) value and the corrected BE(ONIOM2).

### A.5. ONIOM2 Correction

In Table 7 we summarize the single components (namely, BE(CCSD(T), μ-2H₂O) extrapolated, BE(DFT, μ-2H₂O), BE(DFT, ice), and finally the corrected BE(ONIOM2)) of the ONIOM2 correction done with the extrapolation procedure already described in Section 2.

## Appendix B
### Crystalline Adsorption Geometries

The final optimized structures, together with BEs and structural properties for every molecule in our set, are presented in this appendix (see, e.g., Figure 16). For some notable cases, more than one initial structure has been modeled. All energetic quantities are in kelvin (K), while distances are in Å. For every molecule, the gas-phase optimized geometry and its ESP map are also presented. In the latter, the isovalue for the electron density is set equal to $1 \times 10^{-6}$ au, while the values for the ESP (again in au) vary case by case and thus are reported close to the RGB scale legend.

<figure><img src="image_17.png" /><figcaption>**Figure 16.** I) Representation of the hydrogen (H₂) molecule along with its ESP surface, BSSE-corrected binding energies (BE) with (disp) and without (no disp) dispersion. II) Top view of H₂-(010) P-ice surface interaction (unit cell highlighted in blue). III) Detail of the side view.</figcaption></figure>

(The complete figure set (23 images) is available.)
# Appendix C
## Surface Distortion Energy on ASW

In Table 8 we present the distortion energy contributions calculated on the amorphous slab that will be compared with the same contributions calculated on the (010) P-ice model for some notable cases.

<table><caption>**Table 8**
Comparison of the Distortion Energy Contributions (in Kelvin) Calculated on the (010) P-ice Slab and the Mean Value of These Contributions Calculated on the Different Adsorption Cases on the Amorphous Ice Model</caption><thead><tr><th>Species</th><th>δE<sub>s</sub> Crystalline Ice</th><th>⟨δE<sub>s</sub>⟩ Amorphous Ice</th></tr></thead><tbody><tr><td>CH<sub>3</sub>OH</td><td>662</td><td>2790</td></tr><tr><td>H<sub>2</sub>O</td><td>565</td><td>3211</td></tr><tr><td>HCN</td><td>589</td><td>1299</td></tr><tr><td>HCOOH</td><td>878</td><td>1660</td></tr><tr><td>NH<sub>3</sub></td><td>734</td><td>1082</td></tr></tbody></table>

# Appendix D
## Basis Set A-VTZ*

In Tables 9–14 we provide the basis sets in CRYSTAL17 format, used for H, C, N, O, S, and Cl atoms, respectively.

<table><caption>**Table 9**
Basis Set for the H Atom in CRYSTAL17 Format Used in This Work</caption><tbody><tr><td>H</td><td>0</td><td></td><td></td></tr><tr><td>S</td><td>3</td><td>1</td><td></td></tr><tr><td></td><td></td><td>34.06134</td><td>6.03E-03</td></tr><tr><td></td><td></td><td>5.123575</td><td>4.50E-02</td></tr><tr><td></td><td></td><td>1.164663</td><td>0.201897</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.32723</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.103072</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.8</td><td>1</td></tr></tbody></table>

<table><caption>**Table 10**
Basis Set for the C Atom in CRYSTAL17 Format Used in This Work</caption><tbody><tr><td>C</td><td>0</td><td></td><td></td></tr><tr><td>S</td><td>5</td><td>1</td><td></td></tr><tr><td></td><td></td><td>8506.038</td><td>5.34E-04</td></tr><tr><td></td><td></td><td>1275.733</td><td>4.13E-03</td></tr><tr><td></td><td></td><td>290.3119</td><td>2.12E-02</td></tr><tr><td></td><td></td><td>82.0562</td><td>8.24E-02</td></tr><tr><td></td><td></td><td>26.47964</td><td>0.240129</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>9.241459</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>3.364353</td><td>1</td></tr></tbody></table>

<table><caption>**Table 10**
(Continued)</caption><tbody><tr><td>C</td><td>0</td><td></td><td></td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td>S</td><td>1</td><td>0.871742</td><td>1</td></tr><tr><td>S</td><td>1</td><td>0.363524</td><td>1</td></tr><tr><td>S</td><td>1</td><td>0.128731</td><td>1</td></tr><tr><td>P</td><td>4</td><td></td><td></td></tr><tr><td></td><td></td><td>34.7095</td><td>5.33E-03</td></tr><tr><td></td><td></td><td>7.959088</td><td>3.59E-02</td></tr><tr><td></td><td></td><td>2.378697</td><td>0.142003</td></tr><tr><td></td><td></td><td>0.815401</td><td>0.342031</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.289538</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.100848</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>1.6</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.4</td><td>1</td></tr></tbody></table>

<table><caption>**Table 11**
Basis Set for the N Atom in CRYSTAL17 Format Used in This Work</caption><tbody><tr><td>N</td><td>0</td><td></td><td></td></tr><tr><td>S</td><td>5</td><td>1</td><td></td></tr><tr><td></td><td></td><td>11913.42</td><td>-5.23E-04</td></tr><tr><td></td><td></td><td>1786.721</td><td>-4.04E-03</td></tr><tr><td></td><td></td><td>406.5901</td><td>-2.08E-02</td></tr><tr><td></td><td></td><td>114.9253</td><td>-8.12E-02</td></tr><tr><td></td><td></td><td>37.10588</td><td>-0.23871</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>12.97168</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>4.730229</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>1.252518</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.512601</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.179397</td><td>1</td></tr><tr><td>P</td><td>4</td><td></td><td></td></tr><tr><td></td><td></td><td>49.21876</td><td>5.55E-03</td></tr><tr><td></td><td></td><td>11.34894</td><td>3.81E-02</td></tr><tr><td></td><td></td><td>3.428509</td><td>0.149414</td></tr><tr><td></td><td></td><td>1.179951</td><td>0.348982</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.417261</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.142951</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>2</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.5</td><td>1</td></tr></tbody></table>
<table><caption>**Table 12**
Basis Set for the O Atom in CRYSTAL17 Format Used in This Work</caption><thead><tr><th>O</th><th>0</th><th></th><th></th></tr></thead><tbody><tr><td>S</td><td>5</td><td>1</td><td></td></tr><tr><td></td><td></td><td>15902.65</td><td>5.15E-04</td></tr><tr><td></td><td></td><td>2384.954</td><td>3.98E-03</td></tr><tr><td></td><td></td><td>542.7196</td><td>2.05E-02</td></tr><tr><td></td><td></td><td>153.4041</td><td>8.03E-02</td></tr><tr><td></td><td></td><td>49.54572</td><td>0.237668</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>17.33965</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>6.330336</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>1.699588</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.689545</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.23936</td><td>1</td></tr><tr><td>P</td><td>4</td><td></td><td></td></tr><tr><td></td><td></td><td>63.27052</td><td>6.07E-03</td></tr><tr><td></td><td></td><td>14.62331</td><td>4.19E-02</td></tr><tr><td></td><td></td><td>4.448952</td><td>0.161569</td></tr><tr><td></td><td></td><td>1.528151</td><td>0.356828</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.529973</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.175094</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>2.4</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.6</td><td>1</td></tr></tbody></table>

<table><caption>**Table 13**
Basis Set for the S Atom in CRYSTAL17 Format Used in This Work</caption><thead><tr><th>S</th><th>0</th><th></th><th></th></tr></thead><tbody><tr><td>S</td><td>5</td><td></td><td></td></tr><tr><td></td><td></td><td>103954</td><td>2.47E-04</td></tr><tr><td></td><td></td><td>15583.79</td><td>1.92E-03</td></tr><tr><td></td><td></td><td>3546.129</td><td>9.96E-03</td></tr><tr><td></td><td></td><td>1002.681</td><td>4.04E-02</td></tr><tr><td></td><td></td><td>324.9029</td><td>0.130675</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>115.5123</td><td>1</td></tr><tr><td>S</td><td>2</td><td></td><td></td></tr><tr><td></td><td></td><td>44.52821</td><td>0.504036</td></tr><tr><td></td><td></td><td>18.39789</td><td>0.230681</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>5.510068</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>2.125987</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.436919</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.157309</td><td>1</td></tr><tr><td>P</td><td>5</td><td></td><td></td></tr><tr><td></td><td></td><td>606.6937</td><td>2.32E-03</td></tr><tr><td></td><td></td><td>143.507</td><td>1.86E-02</td></tr><tr><td></td><td></td><td>45.74616</td><td>8.60E-02</td></tr><tr><td></td><td></td><td>16.87291</td><td>0.252484</td></tr><tr><td></td><td></td><td>6.63992</td><td>0.446327</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>2.672714</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>1.000009</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.354389</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.116713</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>1.1</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.275</td><td>1</td></tr></tbody></table>

<table><caption>**Table 14**
Basis Set for the Cl Atom in CRYSTAL17 Format Used in This Work</caption><thead><tr><th>Cl</th><th>0</th><th></th><th></th></tr></thead><tbody><tr><td>S</td><td>5</td><td></td><td></td></tr><tr><td></td><td></td><td>117805.8</td><td>2.42E-04</td></tr><tr><td></td><td></td><td>17660.27</td><td>1.87E-03</td></tr><tr><td></td><td></td><td>4018.597</td><td>9.74E-03</td></tr><tr><td></td><td></td><td>1136.223</td><td>3.95E-02</td></tr><tr><td></td><td></td><td>368.1206</td><td>0.127972</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>130.8615</td><td>1</td></tr><tr><td>S</td><td>2</td><td></td><td></td></tr><tr><td></td><td></td><td>50.47901</td><td>0.428741</td></tr><tr><td></td><td></td><td>20.91681</td><td>0.196685</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>6.353139</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>2.494801</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.543359</td><td>1</td></tr><tr><td>S</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.194344</td><td>1</td></tr><tr><td>P</td><td>5</td><td></td><td></td></tr><tr><td></td><td></td><td>681.0688</td><td>2.37E-03</td></tr><tr><td></td><td></td><td>161.1136</td><td>1.89E-02</td></tr><tr><td></td><td></td><td>51.38664</td><td>8.78E-02</td></tr><tr><td></td><td></td><td>18.95851</td><td>0.257074</td></tr><tr><td></td><td></td><td>3.003516</td><td>0.371524</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>7.456529</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>1.060936</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.39452</td><td>1</td></tr><tr><td>P</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.133233</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>1.3</td><td>1</td></tr><tr><td>D</td><td>1</td><td></td><td></td></tr><tr><td></td><td></td><td>0.325</td><td>1</td></tr></tbody></table>

ORCID iDs

Stefano Ferrero https://orcid.org/0000-0001-7819-7657

Lorenzo Zamirri https://orcid.org/0000-0003-0219-6150

Cecilia Ceccarelli https://orcid.org/0000-0001-9664-6292

Arezu Witzel https://orcid.org/0000-0003-0518-944X

Albert Rimola https://orcid.org/0000-0002-9637-4554

Piero Ugliengo https://orcid.org/0000-0001-8886-9832

References

Acharyya, K., & Herbst, E. 2017, ApJ, 850, 105

Al-Halabi, A., & Van Dishoeck, E. 2007, MNRAS, 382, 1648

Åsgeirsson, V., Jónsson, H., & Wikfeldt, K. 2017, JPCC, 121, 1648

Bahr, S., Toubin, C., & Kempter, V. 2008, JChPh, 128, 134712

Balucani, N., Ceccarelli, C., & Taquet, V. 2015, MNRAS: Letters, 449, L16

Bartlett, R. J., & Musial, M. 2007, RvMP, 79, 291

Becke, A. D. 1988, JChPh, 88, 2547

Becke, A. D. 1993, JChPh, 98, 1372

Bennett, C. J., & Kaiser, R. I. 2007, ApJ, 661, 899

Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJ, 570, L101

Boogert, A. A., Gerakines, P. A., & Whittet, D. C. 2015, ARA&A, 53, 541

Borden, W. T., Hoffmann, R., Stuyver, T., & Chen, B. 2017, JACS, 139, 9010

Boys, S. F., & Bernardi, F. 1970, MolPh, 19, 553

Broyden, C. G. 1970, JApMa, 6, 76

Casassa, S. 2000, ChPhL, 321, 1

Casassa, S., Ugliengo, P., & Pisani, C. 1997, JChPh, 106, 8030

Chaabouni, H., Diana, S., Nguyen, T., & Dulieu, F. 2018, A&A, 612, A47

Codella, C., Ceccarelli, C., Bottinelli, S., et al. 2011, ApJ, 744, 164
Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133

Cramer, C. 2002, Essentials of Computational Chemistry (New York: Wiley)

Cuppen, H., Walsh, C., Lamberts, T., et al. 2017, SSRv, 212, 1

Cutini, M., Civalleri, B., Corno, M., et al. 2016, JCTC, 12, 3340

Cutini, M., Civalleri, B., & Ugliengo, P. 2019, ACS Omega, 4, 1838

Cutini, M., Corno, M., & Ugliengo, P. 2017, JCTC, 13, 370

Dapprich, S., Komáromi, I., Byun, K. S., Morokuma, K., & Frisch, M. J. 1999, JMoSt, 461, 1

Das, A., Sil, M., Gorai, P., Chakrabarti, S. K., & Loison, J.-C. 2018, ApJS, 237, 9

Davidson, E. R., & Feller, D. 1986, ChRv, 86, 681

Dohnalek, Z., Kimmel, G. A., Joyce, S. A., et al. 2001, JPhCB, 105, 3747

Douglas, A., & Herzberg, G. 1942, CJRes, 20, 71

Dovesi, R., Erba, A., Orlando, R., et al. 2018, Wiley Interdiscip Rev Comput Mol Sci, 8, e1360

Dulieu, F., Amiaud, L., Congiu, E., et al. 2010, A&A, 512, A30

Dulieu, F., Congiu, E., Noble, J., et al. 2013, NatSR, 3, 1338

Dunning, T. H., Jr 1989, JChPh, 90, 1007

Enrique-Romero, J., Álvarez-Barcia, S., Kolb, F., et al. 2020, MNRAS, 493, 2523

Enrique-Romero, J., Rimola, A., Ceccarelli, C., et al. 2019, ACS Earth Space Chem, 3, 2158

Fayolle, E. C., Balfe, J., Loomis, R., et al. 2016, ApJL, 816, L28

Fletcher, R. 1970, ComJ, 13, 317

Frisch, M., Trucks, G., Schlegel, H., et al. 2009, Gaussian 09 (Wallingford, CT: Gaussian, Inc.), 121, 150

Goldfarb, D. 1970, MaCom, 24, 23

Grimme, S., Antony, J., Ehrlich, S., & Krieg, H. 2010, JChPh, 132, 154104

Grimme, S., Ehrlich, S., & Goerigk, L. 2011, JCoCh, 32, 1456

He, J., Acharyya, K., & Vidali, G. 2016, ApJ, 825, 89

Hiraoka, K., Miyagoshi, T., Takayama, T., Yamamoto, K., & Kihara, Y. 1998, ApJ, 498, 710

Hohenberg, P., & Kohn, W. 1964, PhRv, 136, B864

Hollenbach, D., & Salpeter, E. 1971, ApJ, 163, 155

Horn, A. B., Chesters, M. A., S, M. M. R., & Sodeau, J. R. 1992, Journal of Chemical Society Faraday Transaction, 88, 1077

Jansen, H., & Ros, P. 1969, ChPhL, 3, 140

Jones, A. 2013, A&A, 555, A39

Jones, A., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46

Kama, M., Caux, E., López-Sepulcre, A., et al. 2015, A&A, 574, A107

Karsssemeijer, L., & Cuppen, H. 2014, A&A, 569, A107

Karsssemeijer, L., de Wijs, G., & Cuppen, H. 2014, PCCP, 16, 15630

King, D. A. 1975, SurSc, 47, 384

Kraus, P., & Frank, I. 2018, JPCA, 122, 4894

Lee, C., Yang, W., & Parr, R. G. 1988, PhRvB, 37, 785

Liu, B., & McLean, A. 1973, JChPh, 59, 4557

Martín-Doménech, R., Caro, G. M., Bueno, J., & Goesmann, F. 2014, A&A, 564, A8

McElroy, D., Walsh, C., Markwick, A., et al. 2013, A&A, 550, A36

McGuire, B. A. 2018, ApJS, 239, 17

McKellar, A. 1940, PASP, 52, 187

Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, A&A, 585, A24

Neufeld, D. A., & Wolfire, M. G. 2009, ApJ, 706, 1594

Noble, J., Congiu, E., Dulieu, F., & Fraser, H. 2012, MNRAS, 421, 768

Novoa, J., & Sosa, C. 1995, JPCh, 99, 15837

Oba, Y., Miyauchi, N., Hidaka, H., et al. 2009, ApJ, 701, 464

Oba, Y., Watanabe, N., Hama, T., et al. 2012, ApJ, 749, 67

Öberg, K. I., van Dishoeck, E. F., & Linnartz, H. 2009, A&A, 496, 281

Olanrewaju, B. O., Herring-Captain, J., Grieves, G. A., Aleksandrov, A., & Orlando, T. M. 2011, JPCA, 115, 5936

Pack, J. D., & Monkhorst, H. J. 1977, PhRvB, 16, 1748

Pagani, L., Bourgoin, A., & Lique, F. 2012, A&A, 548, L4

Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2020, ApJ, 897, 56

Papajak, E., Zheng, J., Xu, X., Leverentz, H. R., & Truhlar, D. G. 2011, JCTC, 7, 3027

Pascale, F., Zicovich-Wilson, C. M., Lopez Gejo, F., et al. 2004, JCoCh, 25, 888

Peng, R., Yoshida, H., Chamberlin, R. A., et al. 2010, ApJ, 723, 218

Penteado, E., Walsh, C., & Cuppen, H. 2017, ApJ, 844, 71

Pople, J. A., Gill, P. M., & Handy, N. C. 1995, IJQC, 56, 303

Raut, U., Famá, M., Teolis, B., & Baragiola, R. 2007, JChPh, 127, 204713

Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2019, A&A, 629, A15

Redhead, P. A. 1962, Vacuu, 12, 203

Rimola, A., Civalleri, B., & Ugliengo, P. 2008, Langm, 24, 14027

Rimola, A., Skouteris, D., Balucani, N., et al. 2018, ACS Earth Space Chem, 2, 720

Rimola, A., Taquet, V., Ugliengo, P., Balucani, N., & Ceccarelli, C. 2014, A&A, 572, A70

Schäfer, A., Horn, H., & Ahlrichs, R. 1992, JChPh, 97, 2571

Schilke, P., Phillips, T., & Wang, N. 1995, ApJ, 441, 334

Senevirathne, B., Andersson, S., Dulieu, F., & Nyman, G. 2017, MolAs, 6, 59

Shanno, D. F. 1970, MaCom, 24, 647

Shimonishi, T., Nakatani, N., Furuya, K., & Hama, T. 2018, ApJ, 855, 27

Skouteris, D., Balucani, N., Ceccarelli, C., et al. 2018, ApJ, 854, 135

Smith, R. S., May, R. A., & Kay, B. D. 2016, JPCB, 120, 1979

Song, L., & Kästner, J. 2017, ApJ, 850, 118

Sorrell, W. H. 2001, ApJL, 555, L129

Sure, R., & Grimme, S. 2013, JCoCh, 34, 1672

Svanberg, M., Pettersson, J. B. C., & K. B. 2000, JPCA, 104, 5787

Swings, P., & Rosenfeld, L. 1937, ApJ, 86, 483

Tafalla, M., Myers, P. C., Caselli, P., & Walmsley, C. M. 2004, A&A, 416, 191

Tasker, P. 1979, JPhC, 12, 4977

Tatewaki, H., & Huzinaga, S. 1980, JCoCh, 1, 205

Tentscher, P. R., & Arey, J. S. 2013, JCTC, 9, 1568

Terada, H., & Tokunaga, A. T. 2012, ApJ, 753, 19

Tosoni, S., Pascale, F., Ugliengo, P., et al. 2005, MolPh, 103, 2549

Ugliengo, P., & Damin, A. 2002, ChPhL, 366, 683

Wakelam, V., Loison, J.-C., Herbst, E., et al. 2015, ApJS, 217, 20

Wakelam, V., Loison, J.-C., Mereau, R., & Ruaud, M. 2017, MolAs, 6, 22

Ward, M. D., Hogg, I. A., & Price, S. D. 2012, MNRAS, 425, 1264

Watanabe, N., & Kouchi, A. 2002, ApJL, 571, L173

Watanabe, N., & Kouchi, A. 2008, PrSS, 83, 439

Zamirri, L., Casassa, S., Rimola, A., et al. 2018, MNRAS, 480, 1427

Zamirri, L., Corno, M., Rimola, A., & Ugliengo, P. 2017, ACS Earth Space Chem, 1, 384

Zamirri, L., Pantaleone, S., & Ugliengo, P. 2019a, JChPh, 150, 064702

Zamirri, L., Ugliengo, P., Ceccarelli, C., & Rimola, A. 2019b, ACS Earth Space Chem, 3, 1499

Zhao, Y., & Truhlar, D. G. 2008, Theor Chem Acc, 120, 215

Zicovich-Wilson, C., Pascale, F., Roetti, C., et al. 2004, JCoCh, 25, 1873