# Binding energies of interstellar molecules on crystalline and amorphous
models of water ice by ab-initio calculations

Stefano Ferrero Departament de Química, Universitat Autonoma de Barcelona,
08193 Bellaterra, Catalonia, Spain  Dipartimento di Chimica, Università degli
Studi di Torino, via P. Giuria 7, 10125, Torino, Italy  Lorenzo Zamirri
Dipartimento di Chimica, Università degli Studi di Torino, via P. Giuria 7,
10125, Torino, Italy  Nanostructured Interfaces and Surfaces (NIS) Centre,
Università degli Studi di Torino, via P. Giuria 7, 10125, Torino, Italy
Cecilia Ceccarelli Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
Arezu Witzel Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France  Albert
Rimola Departament de Química, Universitat Autonoma de Barcelona, 08193
Bellaterra, Catalonia, Spain  Piero Ugliengo Dipartimento di Chimica,
Università degli Studi di Torino, via P. Giuria 7, 10125, Torino, Italy

###### Abstract

In the denser and colder ($\leq$20 K) regions of the interstellar medium
(ISM), near-infrared observations have revealed the presence of sub-micron
sized dust grains covered by several layers of H2O-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.

icy dust grains – complex organic molecules – dense molecular clouds –
interstellar medium – ab initio – DFT – surface modeling

††journal: ApJ

tablenum SIXtablenum iint TXFiint

##  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 ($\leq$20–90 K) and densest ($\geq 10^{3}$ cm-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
sub-micron 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 mantles 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 H2O, which is
formed by the hydrogenation of O, O2 and O3 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
$\mu$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, CO2, NH3, CH4, CH3OH and H2CO 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 interstellar ice sample analogous, 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. Oba et al.
2009; Watanabe & Kouchi 2008; 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 which cannot
proceed in the gas phase, such as the formation of H2, the most abundant
molecule in ISM (Hollenbach & Salpeter 1971). Ice surfaces can catalyze
reactions either by behaving as: i) passive third body, this way absorbing
part of the excess of energy released in the surface processes (adsorption
and/or chemical reaction) (e.g. Pantaleone et al. 2020); ii) 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)); and iii) 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 H2CO and CH3OH (e.g.
Watanabe & Kouchi 2002; Rimola et al. 2014; Zamirri et al. 2019b). All three
processes 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, by 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 B _BEs_ for radical
species as they are very reactive. In 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 five hundred 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 which,
in some situations, can overcome the experimental limitations. Many
computational works have so far focused on a few important astrochemical
species like H, H2, N, O, CO, CO2, 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 _BE_ 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 _BE_ 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, whose 4 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 a 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 to
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 modelled 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 fields methods for the
solution of the electronic Schrödinger equation, fully exploiting, if present,
the crystalline 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 non-periodic (molecules) systems with the same
level of accuracy. Furthermore, the definition of the surfaces through the
slab model allows to avoid the 3D fake replica of the slab as forced when
adopting plane waves 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, 14. The needed
density functional integration are 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 by 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 4 _k_
points in the first Brillouin zone. The choice of the numerical values we
assigned to these three computational parameters is fully justified in the
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\text{\times}{10}^{-7}$ Hartree; and maximum components and root-mean-square
of the components of the gradients and atomic displacements vectors,
$4.5\text{\times}{10}^{-4}$ Hartree Bohr-1 and $3\text{\times}{10}^{-4}$
Hartree Bohr-1, and $1.8\text{\times}{10}^{-3}$ Bohr and
$1.2\text{\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 (Becke 1993; Lee et al. 1988),
which has been shown to provide a good level of accuracy for the interaction
energies of non-covalent 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 non-covalent 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 (Kraus & Frank 2018;
Tentscher & Arey 2013). 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 the basis set
(Jansen & Ros 1969; Liu & McLean 1973); ii) the dispersive interactions; 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
$\Gamma$ 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-points formula (NUMDERIV=2 in the CRYSTAL17
code; Dovesi et al. (2018)), based on two displacements of $\pm 0.003$ Å for
each nuclear cartesian coordinates 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 non-
covalent 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, the ice-surface/adsorbate complex to arrive to the
corresponding correction $\Delta ZPE$, as reported in Appendix A.1. From the
$\Delta ZPE$ we corrected the electronic _BE_ for each adsorbates as:
$BE(0)=BE-\Delta ZPE$ and found a good linear correlation _BE(0)=0.854 BE_ ,
as shown in Figure 13 of the Appendix A.1. While the ”fragment frequency”
strategy is fine for computing the $\Delta ZPE$ of the crystalline ice model
due 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 renders the ice structure rather
sensitive to the adsorbate interaction which causes large structural water
molecules rearrangement. This, in turn, alter significantly the whole set of
normal modes and the numerical value of the $\Delta ZPE$ becomes ill-defined.
Nevertheless, considering that the kind of interactions operative for the
crystal ice are of the same nature of those for the amorphous one, we adopted
the same scaling factor 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 Binding energies calculation and Counterpoise correction

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

$\Delta E^{CP}=\Delta E^{*}+\delta E+\Delta E_{L}-BSSE$ (1)

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

$\Delta E^{CP}=-BE$ (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 § 3.2.2. These refinements have been performed through the ONIOM2
approach (Dapprich et al. 1999), dividing the systems in 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 excitations coupled-cluster method added with a
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:

$\displaystyle BE(ONIOM2)=-BE(\emph{Low},\emph{Real})+\Delta BE$ (3)
$\displaystyle\Delta
BE=BE(\emph{High},\emph{Model})-BE(\emph{Low},\emph{Model})$ (4)

The final BE(ONIOM2) is also corrected by the BSSE following the same scheme
described above. Our choices about the _Model_ and _Real system_ will be
extensively justified in § 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 1_ 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 non-periodic 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 twelve atomic layers, is
stoichiometric and has a null electric dipole moment across the z-axis. This
ensure 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 it 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 2x1 supercell. The electrostatic potential maps
(EPMs, see Fig. 1, panels B and C), clearly reveal positive (blue EPM regions)
and negative (red EPM regions) potentials around the dH and dO sites,
respectively.

Figure 1: (010) slab model of P-ice. A) Side view along the _b_ lattice
vector. B) Top view of the 1x1 unit cell ($|$_a_ $|$ = 4.500 Å and $|$_b_ $|$
= 7.078 Å) superimposed to the electrostatic potential map, EPM. C) Top view
of the 2x1 supercell ($|$_a_ $|$ = 8.980 Å, $|$_b_ $|$ = 7.081 Å) superimposed
to the EPM. The iso-surface value for the electron density where the
electrostatic potential is mapped is set equal to 10-6 au. Colour code: +0.02
au (blue, positive), 0.00 au (green, neutral) and -0.02 au (red, negative).

####  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 non-
trivial 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 simulations (MDs) for 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 re-optimized (at B3LYP-D3/A-VTZ* level only) the whole set of ice clusters
and the three most stable clusters, composed by 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 merge of the three clusters was carried
out by matching the dHs regions of one cluster with the dOs 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 by 60 water molecules in the unit
cell, and was further fully optimized (unit cell size and atomic coordinates)
at the HF-3c level, B3LYP-D3/A-VTZ* and M06-2X/A-VTZ* levels of theory. The
three final structures show little differences in the positions of specific
water molecules and, on the whole, the structures are very similar (Fig. 2).
The computed electric dipole moment across the non-periodic direction (1.2,
0.7 and 0.1 Debyes 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 (Fig. 3). The general
characteristics are very similar for the three models: B3LYP-D3 and M06-2X
giving the closest maps. HF-3c tends to enhance the differences between
positive/negative regions due to overpolarization of the electron density
caused by the minimal basis set. “Top” surfaces show a hydrophobic cavity (the
central greenish region, Fig. 3, absent in the P-ice slab, surrounded by dHs
positive spots. “Bottom” surfaces shows several prominent negative regions
(from five dOs) mixed with less prominent positive potentials (due to four
buried dHs).

Figure 2: Side view of the amorphous slab models. The cell parameter _a_ is
highlighted as a blue line. Electric dipole moments $\mu$$\epsilon$ along the
z direction are shown on the right side.

Figure 3: Colour-coded electrostatic potential energy maps (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). dO and dH sites are also
labelled. The iso-surface 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)..

###  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 Fig. 4. For each molecule/surface complex, geometry
optimizations (unit cell plus all atomic coordinates without constraints) were
performed. Initial 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 B3LYP-D3
level of theory. The energetics of the adsorption processes were then computed
according to Eq. (1).

As it 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., H2, N2
and O2) or multipoles moments of higher order (i.e., CH4, see their EPMs in
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
N2, O2 and CH4 cases, interactions are almost repulsive if dispersive
contributions are not accounted for in the total _BE_ (compare _BE_ _disp_
with _BE_ _no disp_ values of Table 1). Therefore, the adsorption is dictated
by dispersive forces, which counterbalance the repulsive electrostatic
interactions. For the H2 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 H2 quadrupole, respectively (see Fig. 16).

Figure 4: Set of molecular and radical species adopted within this work for
the calculation of binding energy on different ice models. O2 is an open-shell
(spin-triplet) species (Borden et al. 2017).

CO, OCS and CO2 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 B) (Zamirri et al. 2017). Thus, although the two negative poles (C and O
atoms) of the quadrupole, both can interact with the positive dH site, the
interaction involving the C atom is energetically slightly favoured 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.

NH3, H2O, HCl, HCN and H2S are all amphiprotic molecules that can both serve
as 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 higher than 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 CH3OH, CH3CN and the three carbonyl-containing
compounds, i.e., H2CO, HCONH2 and HCOOH, all characterized by large molecular
sizes, we adopted the 2x1 supercell (shown in Fig. 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 B). 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• NH2•
CH3• HCO•), since they are of high interest due to their role in the formation
of interstellar compounds (Sorrell 2001; Bennett & Kaiser 2007). OH• and NH2•
form strong H-bonds with the dH and dO sites of the slab, at variance with
CH3• 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_. 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 the Computational Details section, 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 N2, O2, H2O, CH4, CH3CN and CH3• radical. Here, the
_Real system_ is the periodic P-ice slab model without adsorbed species.
Therefore, the BE(_Low_ ,_Real_) term in Eq. (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 by the adsorbed molecule plus $n$ ($n=2,6$; the latter only for the
H2 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 Eq. (3) for any molecular species $\mu$, as BE(CCSD(T), $\mu$-$n$H2O) and
BE(DFT, $\mu$-$n$H2O), respectively.

As CCSD(T) is a wavefunction-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 Jr 1989), here named as cc-pV _N_ Z, 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), $\mu$-$n$H2O) values for _N_ $\to\infty$. Figure 5 shows, using
NH3 as illustrative example, the plot of the BE(CCSD(T), $\mu$-$n$H2O) values
as a function of $1/L^{3}$ 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).

Figure 5: BE(X, $\mu$-$n$H2O) extrapolated value at infinite basis set for the
case of NH3. The dashed-dot blue line represents the _BE_ computed for the
_BE_(DFT, $\mu$-$n$H2O) at DFT//A-VTZ* level (4390 K). Solid red line
represents the linear fit of the BE(CCSD(T), $\mu$-$n$H2O) values (red
squares) calculated with DZ, TZ, QZ and 5Z basis sets. The extrapolated BE(X,
$\mu$-$n$H2O) at infinite basis set is highlighted in red in the fitting
equation (4089 K).

The procedure gives for the extrapolated BE(CCSD(T), $\mu$-$n$H2O) 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 Fig. 5). Very similar
agreement was computed for all considered species as shown in Fig. 6, in which
a very good linear correlation is seen between BE(ONIOM2) and BE(DFT).
Therefore, we can confidently assume the periodic B3LYP-D3/A-VTZ* (closed
shell molecules) or the M06-2X/A-VTZ* (radical species) plain _BE_ values as
reliable and accurate enough and are those actually used in this work.

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:H2; 2:CO; 3:CO2; 4:HCO•;
5:OCS; 6:H2S; 7:HCN; 8:NH2•; 9:H2CO; 10:HCl; 11:OH•; 12:NH3; 13:CH3OH;
14:HCOOH; 15:HCONH2.

####  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 computational costly
when moving from 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 § 2).

To this end, we adopted a composite procedure which 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 discussed for the crystalline ice. We re-optimize each
structure at HF-3c level to check the changes in the structures resulting from
the more approximated method. Then we run a single point energy calculations
at the DFT/A-VTZ* (B3LYP-D3 and M06-2X) levels to evaluate the final _BE_
values. The results obtained are summarized in Fig. 7, showing a very good
linear correlation between the _BE_ values computed as described.

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 empty triangles stand for open-shell and closed-shell
species, respectively.

The largest percentage differences are found for the smallest BEs, that is,
those dominated by dispersion interactions or very weak quadrupolar
interactions (i.e., N2, O2, H2 and CH4) 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 which would have been very
expensive at the full DFT/A-VTZ* level.

###  3.3 BEs on amorphous solid water (ASW)

On the ASW model, due to the presence of different binding sites, a single
_BE_ value is not representative of the whole adsorption processes as it is
the case for almost all adsorbates on the crystalline surface. Therefore, we
computed the _BE_ with the composite DFT//HF-3c procedure (see § 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
setup by hand, following the maximum electrostatic complementarity between the
EPMs (see Fig. 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 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 10091 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 a twice depending on the site and that the
value on the crystalline ice can also be substantially different to that on
the ASW.

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. 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_.

Species | (010) P-ice crystalline slab DFT//DFT |  | (010) P-ice crystalline slab DFT//HF-3c
---|---|---|---
_BE disp_ | _BE no disp_ | -disp(%) |  | _BE disp_ | _BE no disp_ | -disp(%)
H2 | 1191 | 565 | 625(53) |  | 926 | 241 | 686(74)
O2 | 1022 | -373 | 1034(137) |  | 794 | -84 | 878(110)
N2 | 1564 | -72 | 1636(104) |  | 1455 | -180 | 1636(160)
CH4 | 1684 | -229 | 1912(113) |  | 1912 | -349 | 2261(118)
CO | 2357 | 698 | 1660(71) |  | 1948 | 60 | 1888 (97)
CO2 | 3440 | 1540 | 1900(55) |  | 3007 | 938 | 2069(69)
OCS | 3476 | 120 | 3356(97) |  | 3187 | 265 | 2923(92)
HCl | 6507 | 4402 | 2093(32) |  | 6314 | 3488 | 2237(39)
HCN | 5124 | 3067 | 2057(29) |  | 5725 | 3271 | 3043(48)
H2O | 8431 | 6844 | 1588(19) |  | 8431 | 6808 | 1612(19)
H2S | 5677 | 3380 | 2297(40) |  | 5232 | 3199 | 2105(40)
NH3 | 7373 | 5533 | 1852(25) |  | 7301 | 5484 | 1816(25)
CH3CN | 7553 | 4450 | 3103(41) |  | 6916 | 3259 | 2598(44)
CH3OH | 8684 | 6014 | 2670(31) |  | 8648 | 6026 | 2237(27)
H2CO-SC1 | 5869 | 3885 | 1985(34) |  | 5773 | 4053 | 2369(37)
H2CO-SC2 | 6375 | 3692 | 2682(42) |  | 6423 | 3716 | 2057(36)
HCONH2-SC1 | 9610 | 6459 | 3151(33) |  | 9321 | 6158 | 3163(34)
HCONH2-SC2 | 10079 | 6483 | 3608(36) |  | 9634 | 6074 | 3560(37)
HCOOH | 9526 | 7325 | 2189(23) |  | 9297 | 7168 | 2117(23)
HCOOH-SC | 9442 | 7301 | 2021(21) |  | 9405 | 7541 | 1864(20)
OH• | 6543* |  |  |  | 6795* |  |
HCO• | 3476* |  |  |  | 3548* |  |
CH3• | 2562* |  |  |  | 2598* |  |
NH2• | 6038* |  |  |  | 6050* |  |

  * •

*Notes: For radical species (energy at M06-2X level) we cannot discern between _disp_ and _no disp_ data.

Figure 8: Comparison of the final optimized geometries for CH3OH and HCONH2
(as illustrative examples) on the crystalline ice (§ 3.1.1) and on the ASW (§
3.1.2). The _BE_ values (in Kelvin) are reported in each plot.

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 H2, O2,
N2, CH4, CO, CO2, 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 becomes dissociated at the amorphous one. Further
details about the case of HCl are reported in section A.2 of the Appendix.
Probably this will not occur for HF, not considered here, which is expected to
be molecularly adsorbed on both ices due 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.

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).

Some adsorbates show similar trends in the BEs, despite their different
chemical nature. This is shown in Fig. 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 H2CO and HCOOH are very similar (in their relative values), and those for
CH3OH and HCONH2 show some similarities, despite the large difference in the
chemical functionality.

Table 2: _BE_ values ($\mathrm{K}$) calculated with DFT//HF-3c method for every case on the amorphous slab model, where the ZPE correction has not been added. | Amorphous ice _BE_ values
---|---
Species | Case 1 | Case 2 | Case 3 | Case 4 | Case 5 | Case 6 | Case 7 | Case 8
H2 | 469 | 505 | 277 | 361 | 265 |  |  |
O2 | 818 | 854 | 854 | 529 | 337 |  |  |
N2 | 1347 | 1708 | 1311 | 1191 | 890 |  |  |
CH4 | 1323 | 1960 | 1636 | 1467 | 1070 |  |  |
CO | 1816 | 2189 | 1540 | 1527 | 1299 |  |  |
CO2 | 2863 | 3452 | 2538 | 2550 | 1744 |  |  |
OCS | 3404 | 2971 | 2670 | 2562 | 1780 |  |  |
HCN | 2923 | 5124 | 3620 | 5136 | 3271 | 5857 | 6146 | 7421
H2O | 7156 | 5845 | 6014 | 5689 | 4222 |  |  |
H2S | 2814 | 3909 | 3151 | 3560 | 2682 |  |  |
NH3 | 8840 | 5268 | 6820 | 5930 | 6579 | 6098 | 5052 |
CH3CN | 8960 | 5857 | 5617 | 5557 | 6483 |  |  |
CH3OH | 6531 | 4414 | 6519 | 6362 | 5208 | 6663 | 5509 | 10091
H2CO | 3596 | 4258 | 4174 | 4775 | 5268 | 4594 | 3873 | 7253
HCONH2 | 12833 | 7481 | 10344 | 6820 | 8467 | 7072 | 6783 | 7313
HCOOH | 7577 | 7409 | 8515 | 12364 | 6302 | 7204 | 6639 | 11354
OH• | 6230 | 1816 | 4955 | 2754 | 5076 |  |  |
HCO• | 2694 | 2057 | 1540 | 3608 | 1672 |  |  |
CH3• | 1708 | 1936 | 1299 |  |  |  |  |
NH2• | 4354 | 3716 | 4402 | 3368 |  |  |  |
Table 3: Summary of our computed _BEs_ and comparison with data from the
literature. The first column reports the species, the columns 2 to 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, column 7 and 8 the values
in the two astrochemical databases KIDA and UMIST (see text), and the last
columns the values measured in different experiments. Units are in K and the
references are listed in the notes below.

| BEs from this work |  | BEs from literature
---|---|---|---
| Crystalline ice | ASW |  | Computed | Databases | Experiments
Species | BE(0) _disp_ | Min | Max |  | Das(a) | Wakelam(b) | UMIST(c) | KIDA(d) | Penteado(e) | Others
H2 | 790 | 226 | 431 |  | 545 | 800 | 430 | 440 | 480 $\pm$ 10 | 322-505(f)
O2 | 677 | 287 | 729 |  | 1352 | 1000 | 1000 | 1200 | 914-1161 | 920-1520(f)
N2 | 1242 | 760 | 1458 |  | 1161 | 1100 | 790 | 1100 | 1200 | 790-1320(f1)
|  |  |  |  |  |  |  |  |  | 900-1800(f1)
CH4 | 1633 | 914 | 1674 |  | 2321 | 800 | 1090 | 960 | 1370 | 960-1947(f,g,h)
CO | 1663 | 1109 | 1869 |  | 1292 | 1300 | 1150 | 1300 | 863-1420 | 870-1600(f)
|  |  |  |  |  |  |  |  |  | 980-1940(f1)
CO2 | 2568 | 1489 | 2948 |  | 2352 | 3100 | 2990 | 2600 | 2236-2346 |
OCS | 2722 | 1520 | 2907 |  | 1808 | 2100 | 2888 | 2400 | 2325** | 2430(i)
HCl | 5557 | (l) | (l) |  | 4104 | 4800 |  | 5172 |  | 5172(m)
HCN | 6392 | 2496 | 6337 |  | 2352 | 3500 | 2050 | 3700 |  |
H2O | 7200 | 3605 | 6111 |  | 4166 | 4600 | 4800 | 5600 | 4815-5930 |
H2S | 4468 | 2291 | 3338 |  | 3232 | 2500-2900 | 2743 | 2700 | 2296** |
NH3 | 6235 | 4314 | 7549 |  | 5163 | 5600 | 5534 | 5500 | 2715** |
CH3CN | 5906 | 4745 | 7652 |  | 3786 | 4300 | 4680 | 4680 | 3790** |
CH3OH | 7385 | 3770 | 8618 |  | 4511 | 4500-5100 | 4930 | 5000 | 3820** | 3700-5410(n,o,p)
H2CO | 5187 | 3071 | 6194 |  | 3242 | 5100 | 2050 | 4500 | 3260$\pm$60 |
HCONH2 | 8104 | 5793 | 10960 |  |  | 6300 | 5556 |  |  | 7460-9380(q)
HCOOH | 7991 | 5382 | 10559 |  | 3483 |  | 5000 | 5570 | 4532** |
OH• | 5588 | 1551 | 5321 |  | 3183 | 3300-5300 | 2850 | 4600 | 1656-4760 |
HCO• | 2968 | 1315 | 3081 |  | 1857 | 2300-2700 | 1600 | 2400 |  |
CH3• | 2188 | 1109 | 1654 |  | 1322 | 2500 | 1175 | 1600 |  |
NH2• | 5156 | 2876 | 4459 |  | 3240 | 2800-4500 | 3956 | 3200 |  |

  * •

*Notes: (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. (2011b) (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; ** Results estimated from the work of Collings et al. (2004), reported in Table 2 of Penteado et al. (2017)

##  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
two (Table 3). This has already 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
_BEs(0)_ with those in the literature § 4.1 and for the astrophysical
implications § 4.2. We will discuss these two aspects separately in the next
two sections.

Figure 10: Spider graphs of the DFT//HF-3c _BE_ values (in Kelvin) calculated
on the same adsorption eight sites of the ASW for the H2CO (red), CH3OH
(black), HCONH2 (green) and HCOOH (blue) molecules. The _BE_ values scale goes
from 0 Kelvin (center of the graph) to 14433 Kelvin (vertices of the polygon)
in steps of 4811 Kelvin. Labeling of dH and dO sites is referring to Fig. 3

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., Wakelam et al.
and reported in the KIDA and UMIST databases (McElroy et al. 2013; Wakelam et
al. 2015, 2017; Das et al. 2018).

###  4.1 Comparison _BE_ values in 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 binding energy (_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\longrightarrow 0$). The correct comparison with experiments would
therefore be with microcalorimetric measurements at the zero-limit adsorbate
coverage. In astrochemical laboratories, Temperature Programmed Desorption
(TPD) is, instead, the method of choice and is related to the desorption
activation energy (_DAE_). _DAE_ derives indirectly from the TPD peaks through
the Readhead’s method (Redhead 1962), or more sophisticated techniques. TPD
usually start 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 (Collings et al. 2004). Therefore, _BE_ reported in the review by Penteado
et al (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_{2}O)]BE(H_{2}O)$, in which $T_{des}(X)$ is the
desorption temperature of the $X$ species contrasted with that of water
$T_{des}(H_{2}O)$ to arrive to 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
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 N2 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
_BE_ on amorphous ice models are in reasonable agreement with that measured by
He et al. for CO (1109-1869 K versus 1940 K), and on the lower side for N2
(760-1455 K versus 1800 K). It is worth-noting that the comparison is much
better when referring to the non-porous amorphous ice measurements by the same
authors: for CO and N2 the measured _BE_ values are 1600 K and 1320 K, indeed
well bracketed by our _BE(0)_.

Data from Penteado et al (Penteado et al. 2017) extracted from the TPD of
Collings et al (Collings et al. 2004) for NH3 adsorbed on ice layer gave a
_BE_ of 2715 K. This value is, however, identical to that from the TPD of NH3
adsorbed on gold surface (no water ice), proving that that _BE_ is relative to
the NH3/NH3 lateral interaction within the adsorbed NH3 multilayer and not due
to the interaction with the ice surface. For reasons explained in Collings et
al (Collings et al. 2004) the _BE_ of NH3 on ice is assumed of the same order
of 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 brackets the experimental ones. The computed highest values
emphasize the H-bond acceptor capability of NH3 occurring on few specific
sites characterized by very high electrostatic potential, only important for
very low NH3 coverage, not easily accessible in the TPD experiments. For the
H2O 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 a relatively
different computed and measured BEs: O2. For O2, 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 H2,
our computed lowest _BEs(0)_ are within the range of the measured ones, but we
predict sites where _BEs(0)_ are larger, which may have important
astrophysical implications (§ 4.2). Finally, for H2 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 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 none 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 to fully considering the hydrogen bond cooperativity, which
will enhance the strength of the interaction with adsorbates at the terminal
dangling hydrogen 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 Fig. 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 CH3CN 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 H2
and O2.

####  4.1.3 Comparison with values in astrochemical databases

Two databases list the _BE_ of the species used by the astrochemical models:
the KIDA (Kinetic Database for Astrochemistry,
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 Fig. 11. The general remarks that we wrote for the comparison with
the literature experimental and theoretical values (§ 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 binding energies 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) 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 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 of the ice,
namely exposed to the void, can sublimate (and not the entire bulk). In this
toy model, we considered ten layers for an icy grain mantle, where the bottom
five are made entirely of water and each of the top five layers contains 80%
of water and 20% of methanol. The methanol molecules with different _BEs_ are
distributed randomly on each layer, with a 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 10091 K) and, at each layer, there will be 12.5% of 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% molecule 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 Fig. 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 surprising, the introduction of multiple _BEs_ produces multiple peaks of
desorption, both for water and methanol. Figure 12 shows that the water
desorption rate has a small peak at $\sim$75 K, a larger one at $\sim$120 K,
then another at $\sim$140 K and, finally, a last peak at $\sim$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_.

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 by ten layers: the bottom five layers contain water while the top
fives contain 80% of water and 20% of 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 normalised to one.

####  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 10091 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 towards hot cores and hot corinos, or, finally, towards
prestellar cores.

####  4.2.3 Comments on N2, CO and HCl

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

N2 and CO: Our computations show that the _BE_ of CO is definitively larger
than that of N2, 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 data bases the difference is 200 and 360 K in KIDA and UMIST,
respectively. This difference very likely can explain why observations detect
N2H+, which is formed in the gas phase from N2, 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 which 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 out of the scope of the present work. Here we want just 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 binding energy is about 500 K larger than that of N2, 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 so far carried out 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, this latter obviously unobservable. This
would help solving the mystery of HCl not being observed in gas-phase.
Furthermore, if that would be the case, the population of the chemically
reactive atomic Cl will be increased, with 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 binding energies (_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 models. We relied on density functional theory (DFT)
calculations, using the B3LYP-D3 and M06-2X widely used functionals. This
approach was further validated by an ONIOM-like correction at 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 binding energies computed at the crystalline
ice surface at DFT//DFT and DFT//HF-3c level, which were found 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, both from 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 laboratory, with the
exception of O2 and, at less extent, H2. On the contrary, previous
computations of _BEs_ , which considered a very small ($\leq 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 which
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 _BEs_
values.

We discussed some astrophysical implications, showing that the multiple
computed _BEs_ give rise to a complex process of interstellar iced mantles
desorption, with multiple peaks as a function of the temperature that depends
(also) on the ice structure. Our new computations do not allow to estimate how
the _BEs_ are distributed for each molecule as only a statistical study is
necessary for that. The new (multiple) _BEs_ of N2 and CO might explain why
N2H+ deplete 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 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. SF, LZ and PU 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. AR 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 programme under the
Marie Sklodowska-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 programme, 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.

## References

  * Acharyya & Herbst (2017) Acharyya, K., & Herbst, E. 2017, ApJ, 850, 105
  * Al-Halabi & Van Dishoeck (2007) Al-Halabi, A., & Van Dishoeck, E. 2007, MNRAS, 382, 1648
  * Ásgeirsson et al. (2017) Ásgeirsson, V., Jónsson, H., & Wikfeldt, K. 2017, The Journal of Physical Chemistry C, 121, 1648
  * Bahr et al. (2008) Bahr, S., Toubin, C., & Kempter, V. 2008, The Journal of chemical physics, 128, 134712
  * Balucani et al. (2015) Balucani, N., Ceccarelli, C., & Taquet, V. 2015, MNRAS: Letters, 449, L16
  * Bartlett & Musiał (2007) Bartlett, R. J., & Musiał, M. 2007, Rev Mod Phys, 79, 291
  * Becke (1988) Becke, A. D. 1988, J Chem Phys, 88, 2547
  * Becke (1993) —. 1993, The Journal of Chemical Physics, 98, 1372
  * Bennett & Kaiser (2007) Bennett, C. J., & Kaiser, R. I. 2007, ApJ, 661, 899
  * Bergin et al. (2002) Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJ, 570, L101
  * Boogert et al. (2015) Boogert, A. A., Gerakines, P. A., & Whittet, D. C. 2015, Annu Rev Astron Astrophys, 53, 541
  * Borden et al. (2017) Borden, W. T., Hoffmann, R., Stuyver, T., & Chen, B. 2017, JACS, 139, 9010
  * Boys & Bernardi (1970) Boys, S. F., & Bernardi, F. 1970, Mol Phys, 19, 553
  * Broyden (1970) Broyden, C. G. 1970, IMA J Appl Math, 6, 76
  * Casassa (2000) Casassa, S. 2000, Chem Phys Lett, 321, 1
  * Casassa et al. (1997) Casassa, S., Ugliengo, P., & Pisani, C. 1997, J Chem Phys, 106, 8030
  * Chaabouni et al. (2018) Chaabouni, H., Diana, S., Nguyen, T., & Dulieu, F. 2018, A&A, 612, A47
  * Codella et al. (2011) Codella, C., Ceccarelli, C., Bottinelli, S., et al. 2011, ApJ, 744, 164
  * Collings et al. (2004) Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133
  * Cramer (2002) Cramer, C. 2002, Essentials of computational chemistry (England, Wiley
  * Cuppen et al. (2017) Cuppen, H., Walsh, C., Lamberts, T., et al. 2017, SSRev, 212, 1
  * Cutini et al. (2016) Cutini, M., Civalleri, B., Corno, M., et al. 2016, JCTC, 12, 3340
  * Cutini et al. (2019) Cutini, M., Civalleri, B., & Ugliengo, P. 2019, ACS omega, 4, 1838
  * Cutini et al. (2017) Cutini, M., Corno, M., & Ugliengo, P. 2017, JCTC, 13, 370
  * Dapprich et al. (1999) Dapprich, S., Komáromi, I., Byun, K. S., Morokuma, K., & Frisch, M. J. 1999, J Mol Struc:THEOCHEM, 461, 1
  * Das et al. (2018) Das, A., Sil, M., Gorai, P., Chakrabarti, S. K., & Loison, J.-C. 2018, ApJS, 237, 9
  * Davidson & Feller (1986) Davidson, E. R., & Feller, D. 1986, Chem Rev, 86, 681
  * Dohnalek et al. (2001) Dohnalek, Z., Kimmel, G. A., Joyce, S. A., et al. 2001, J Phys Chem B, 105, 3747
  * Douglas & Herzberg (1942) Douglas, A., & Herzberg, G. 1942, Canadian Journal of Research, 20, 71
  * Dovesi et al. (2018) Dovesi, R., Erba, A., Orlando, R., et al. 2018, Wiley Interdiscip Rev Comput Mol Sci, 8, e1360
  * Dulieu et al. (2010) Dulieu, F., Amiaud, L., Congiu, E., et al. 2010, A&A, 512, A30
  * Dulieu et al. (2013) Dulieu, F., Congiu, E., Noble, J., et al. 2013, Sci Rep, 3, 1338
  * Dunning Jr (1989) Dunning Jr, T. H. 1989, J Chem Phys, 90, 1007
  * Enrique-Romero et al. (2019) Enrique-Romero, J., Rimola, A., Ceccarelli, C., et al. 2019, ACS Earth Space Chem, 3, 2158
  * Enrique-Romero et al. (2020) Enrique-Romero, J., Álvarez-Barcia, S., Kolb, F., et al. 2020, MNRAS, 493, 2523
  * Fayolle et al. (2016) Fayolle, E. C., Balfe, J., Loomis, R., et al. 2016, ApJ, 816, L28
  * Fletcher (1970) Fletcher, R. 1970, Comput J, 13, 317
  * Frisch et al. (2009) Frisch, M., Trucks, G., Schlegel, H., et al. 2009, Inc., Wallingford CT, 121, 150
  * Goldfarb (1970) Goldfarb, D. 1970, Math Comput, 24, 23
  * Grimme et al. (2010) Grimme, S., Antony, J., Ehrlich, S., & Krieg, H. 2010, J Chem Phys, 132, 154104
  * Grimme et al. (2011) Grimme, S., Ehrlich, S., & Goerigk, L. 2011, J Comput Chem, 32, 1456
  * He et al. (2016) He, J., Acharyya, K., & Vidali, G. 2016, ApJ, 825, 89
  * Hiraoka et al. (1998) Hiraoka, K., Miyagoshi, T., Takayama, T., Yamamoto, K., & Kihara, Y. 1998, ApJ, 498, 710
  * Hohenberg & Kohn (1964) Hohenberg, P., & Kohn, W. 1964, Phys Rev, 136, B864
  * Hollenbach & Salpeter (1971) Hollenbach, D., & Salpeter, E. 1971, ApJ, 163, 155
  * Horn et al. (1992) Horn, A. B., Chesters, M. A., S, M. M. R., & Sodeau, J. R. 1992, Journal of Chemical Society Faraday Transaction, 88, 1077
  * Jansen & Ros (1969) Jansen, H., & Ros, P. 1969, Chem Phys Lett, 3, 140
  * Jones (2013) Jones, A. 2013, A&A, 555, A39
  * Jones et al. (2017) Jones, A., Köhler, M., Ysard, N., Bocchio, M., & Verstraete, L. 2017, A&A, 602, A46
  * Kama et al. (2015) Kama, M., Caux, E., López-Sepulcre, A., et al. 2015, A&A, 574, A107
  * Karssemeijer & Cuppen (2014) Karssemeijer, L., & Cuppen, H. 2014, A&A, 569, A107
  * Karssemeijer et al. (2014) Karssemeijer, L., de Wijs, G., & Cuppen, H. 2014, PCCP, 16, 15630
  * King (1975) King, D. A. 1975, Surface Science, 47, 384
  * Kraus & Frank (2018) Kraus, P., & Frank, I. 2018, J Phys Chem A, 122, 4894
  * Lee et al. (1988) Lee, C., Yang, W., & Parr, R. G. 1988, Phys Rev B, 37, 785
  * Liu & McLean (1973) Liu, B., & McLean, A. 1973, J Chem Phys, 59, 4557
  * Martín-Doménech et al. (2014) Martín-Doménech, R., Caro, G. M., Bueno, J., & Goesmann, F. 2014, Astronomy & Astrophysics, 564, A8
  * McElroy et al. (2013) McElroy, D., Walsh, C., Markwick, A., et al. 2013, A&A, 550, A36
  * McGuire (2018) McGuire, B. A. 2018, ApJS, 239, 17
  * McKellar (1940) McKellar, A. 1940, Publications of the Astronomical Society of the Pacific, 52, 187
  * Minissale et al. (2016) Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, Astronomy & Astrophysics, 585, A24
  * Neufeld & Wolfire (2009) Neufeld, D. A., & Wolfire, M. G. 2009, ApJ, 706, 1594
  * Noble et al. (2012) Noble, J., Congiu, E., Dulieu, F., & Fraser, H. 2012, MNRAS, 421, 768
  * Novoa & Sosa (1995) Novoa, J., & Sosa, C. 1995, The Journal of Physical Chemistry, 99, 15837
  * Oba et al. (2009) Oba, Y., Miyauchi, N., Hidaka, H., et al. 2009, ApJ, 701, 464
  * Oba et al. (2012) Oba, Y., Watanabe, N., Hama, T., et al. 2012, ApJ, 749, 67
  * Öberg et al. (2009) Öberg, K. I., van Dishoeck, E. F., & Linnartz, H. 2009, A&A, 496, 281
  * Olanrewaju et al. (2011a) Olanrewaju, B. O., Herring-Captain, J., A, G. G., Aleksandrov, A., & M, O. T. 2011a, The Journal of Physical Chemistry A, 115, 5936
  * Olanrewaju et al. (2011b) Olanrewaju, B. O., Herring-Captain, J., Grieves, G. A., Aleksandrov, A., & Orlando, T. M. 2011b, The Journal of Physical Chemistry A, 115, 5936
  * Pack & Monkhorst (1977) Pack, J. D., & Monkhorst, H. J. 1977, Phys Rev B, 16, 1748
  * Pagani et al. (2012) Pagani, L., Bourgoin, A., & Lique, F. 2012, A&A, 548, L4
  * Pantaleone et al. (2020) Pantaleone, S., Enrique-Romero, J., Ceccarelli, C., et al. 2020, Chemical desorption versus energy dissipation: insights from ab-initio molecular dynamics of HCO formation, , , arXiv:2004.11758
  * Papajak et al. (2011) Papajak, E., Zheng, J., Xu, X., Leverentz, H. R., & Truhlar, D. G. 2011, JCTC, 7, 3027
  * Pascale et al. (2004) Pascale, F., Zicovich-Wilson, C. M., Lopez Gejo, F., et al. 2004, J Comput Chem, 25, 888
  * Peng et al. (2010) Peng, R., Yoshida, H., Chamberlin, R. A., et al. 2010, ApJ, 723, 218
  * Penteado et al. (2017) Penteado, E., Walsh, C., & Cuppen, H. 2017, ApJ, 844, 71
  * Pople et al. (1995) Pople, J. A., Gill, P. M., & Handy, N. C. 1995, Int J Quantum Chem, 56, 303
  * Raut et al. (2007) Raut, U., Famá, M., Teolis, B., & Baragiola, R. 2007, The Journal of chemical physics, 127, 204713
  * Redaelli et al. (2019) Redaelli, E., Bizzocchi, L., Caselli, P., et al. 2019, A&A, 629, A15
  * Redhead (1962) Redhead, P. A. 1962, Vacuum, 12, 203
  * Rimola et al. (2008) Rimola, A., Civalleri, B., & Ugliengo, P. 2008, Langmuir, 24, 14027
  * Rimola et al. (2018) Rimola, A., Skouteris, D., Balucani, N., et al. 2018, ACS Earth Space Chem, 2, 720
  * Rimola et al. (2014) Rimola, A., Taquet, V., Ugliengo, P., Balucani, N., & Ceccarelli, C. 2014, A&A, 572, A70
  * Schäfer et al. (1992) Schäfer, A., Horn, H., & Ahlrichs, R. 1992, J Chem Phys, 97, 2571
  * Schilke et al. (1995) Schilke, P., Phillips, T., & Wang, N. 1995, ApJ, 441, 334
  * Senevirathne et al. (2017) Senevirathne, B., Andersson, S., Dulieu, F., & Nyman, G. 2017, Mol Astrophys, 6, 59
  * Shanno (1970) Shanno, D. F. 1970, Math Comput, 24, 647
  * Shimonishi et al. (2018) Shimonishi, T., Nakatani, N., Furuya, K., & Hama, T. 2018, ApJ, 855, 27
  * Skouteris et al. (2018) Skouteris, D., Balucani, N., Ceccarelli, C., et al. 2018, ApJ, 854, 135
  * Smith et al. (2016) Smith, R. S., May, R. A., & Kay, B. D. 2016, The Journal of Physical Chemistry B, 120, 1979
  * Song & Kästner (2017) Song, L., & Kästner, J. 2017, ApJ, 850, 118
  * Sorrell (2001) Sorrell, W. H. 2001, ApJ, 555, L129
  * Sure & Grimme (2013) Sure, R., & Grimme, S. 2013, J Comput Chem, 34, 1672
  * Svanberg et al. (2000) Svanberg, M., Pettersson, J. B. C., & K, B. 2000, The Journal of Physical Chemistry A, 104, 5787
  * Swings & Rosenfeld (1937) Swings, P., & Rosenfeld, L. 1937, ApJ, 86, 483
  * Tafalla et al. (2004) Tafalla, M., Myers, P. C., Caselli, P., & Walmsley, C. M. 2004, A&A, 416, 191
  * Tasker (1979) Tasker, P. 1979, Journal of Physics C: Solid State Physics, 12, 4977
  * Tatewaki & Huzinaga (1980) Tatewaki, H., & Huzinaga, S. 1980, J Comput Chem, 1, 205
  * Tentscher & Arey (2013) Tentscher, P. R., & Arey, J. S. 2013, JCTC, 9, 1568
  * Terada & Tokunaga (2012) Terada, H., & Tokunaga, A. T. 2012, ApJ, 753, 19
  * Tosoni et al. (2005) Tosoni, S., Pascale, F., Ugliengo, P., et al. 2005, Mol Phys, 103, 2549
  * Ugliengo & Damin (2002) Ugliengo, P., & Damin, A. 2002, Chem Phys Lett, 366, 683
  * Wakelam et al. (2017) Wakelam, V., Loison, J.-C., Mereau, R., & Ruaud, M. 2017, Mol Astrophys, 6, 22
  * Wakelam et al. (2015) Wakelam, V., Loison, J.-C., Herbst, E., et al. 2015, ApJS, 217, 20
  * Ward et al. (2012) Ward, M. D., Hogg, I. A., & Price, S. D. 2012, Monthly Notices of the Royal Astronomical Society, 425, 1264
  * Watanabe & Kouchi (2002) Watanabe, N., & Kouchi, A. 2002, ApJ, 571, L173
  * Watanabe & Kouchi (2008) —. 2008, Prog Surf Sci, 83, 439
  * Zamirri et al. (2018) Zamirri, L., Casassa, S., Rimola, A., et al. 2018, MNRAS, 480, 1427
  * Zamirri et al. (2017) Zamirri, L., Corno, M., Rimola, A., & Ugliengo, P. 2017, ACS Earth Space Chem, 1, 384
  * Zamirri et al. (2019a) Zamirri, L., Pantaleone, S., & Ugliengo, P. 2019a, J Chem Phys, 150, 064702
  * Zamirri et al. (2019b) Zamirri, L., Ugliengo, P., Ceccarelli, C., & Rimola, A. 2019b, ACS Earth Space Chem, 3, 1499
  * Zhao & Truhlar (2008) Zhao, Y., & Truhlar, D. G. 2008, Theor Chem Acc, 120, 215
  * Zicovich-Wilson et al. (2004) Zicovich-Wilson, C., Pascale, F., Roetti, C., et al. 2004, J Comput Chem, 25, 1873

##  Appendix A Computational details

In CRYSTAL17, the multielectron wavefunction is built as a Slater determinant
of crystalline/molecular orbitals which are linear combinations of localized
functions on the different atoms of the structure which are called Atomic
Orbitals (AOs). In a similar manner, the AOs are constructed by linear
combinations of localized Gaussian functions which 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 Binding energies, Counterpoise and zero point energy corrections

In a periodic treatment of surface adsorption phenomena one of the most
relevant energy value, useful to describe the interacting system, is the
binding energy _BEs_ which is related to the interaction energy $\Delta E$, so
that:

$BE=-\Delta E$ (A1)

The binding energy 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)$ (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\textsubscript{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 h

$\displaystyle BE$ $\displaystyle=BE^{*}-\delta E_{S}-\delta E_{M}$ (A3)
$\displaystyle\delta E_{S}$ $\displaystyle=E(S//SM)-E(S//S)$ (A4)
$\displaystyle\delta E_{M}$ $\displaystyle=E(M//SM)-E\textsubscript{m}(M//M)$
(A5) $\displaystyle BE^{*}$ $\displaystyle=[E(S//SM)+E(M//SM)]-E(SM//SM)$ (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 purely molecule’s deformation energy can be computed as:

$\Delta E_{M}=E_{m}(M//SM)-E_{m}(M//M)$ (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)$ (A8)

and can be either positive (repulsion) or negative (attraction). With those
positions, the $BE^{*}$ interaction energy is then deformation and lateral
interactions 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 the same geometry assumed in the $SM$ system. The above _BE_
definition can be easily recast to include the basis set superposition error
(_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 $BE^{C}$ is then:

$\displaystyle BE^{C}$ $\displaystyle=BE^{*C}-\delta E_{S}-\Delta E_{M}-\Delta
E_{L}^{C}$ (A9) $\displaystyle BE^{*C}$
$\displaystyle=[E(S[M]//SM)+E([S]M//SM)]-E(SM//SM)$ (A10) $\displaystyle BSSE$
$\displaystyle=BE-BE^{C}$ (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 $BE^{*C}<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 to the definition of the zero point correct binding energy $BE(0)$
as:

$\displaystyle BE(0)$
$\displaystyle=[E_{m}(M//M)+E(S//S)]-E(SM//SM)+[ZPE(M)+ZPE(S)-ZPE(SM)]$ (A12)
$\displaystyle BE(0)$ $\displaystyle=BE-\Delta ZPE$ (A13) $\displaystyle\Delta
ZPE$ $\displaystyle=ZPE(SM)-ZPE(M)-ZPE(S)$ (A14)

In this work the $\Delta ZPE$ has been evaluated for the crystalline ice cases
only and the same correction was adopted also to correct the $BE$ for the the
amorphous ice model for the reason described in the Computational Details
section.

Figure 13: Linear regression between _BE_ and _BE(0)_ (zero point corrected)
values computed for the considered adsorbates on the crystalline ice

###  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 due 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_ (Svanberg et al. 2000) when enough dangling hydrogen 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 early
found by Horn _et al_ (Horn et al. 1992) through RAIRS spectroscopy of DCl
adsorption on thin D2O film, and by Olanrewaju _et al_ (Olanrewaju et al.
2011a) through thermal and electron-stimulated desorption experiments.

Figure 14: Enlarged views of HCl adsorbed on a crystalline and amorphous ice
models

###  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
(Dovesi et al. 2018) in the code) as described in §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 $\Delta E$ as define by Equation (A1) being the
threshold. For this benchmark calculations, $BSSE$, distortion and lateral
interaction contributions have not been computed. Results are resumed in what
follows.

Table 4: Optimization of the SHRINK parameter. Tolerances of integrals (TOLINTEG) and integration grid (XLGRID) are fixed to default values. SHRINK |  _k_ points in BZ |  $\Delta E$ ($\mathrm{K}$)
---|---|---
2 2 | 4 | 7873
3 3 | 5 | 7889
4 4 | 10 | 7889
5 5 | 13 | 7890
6 6 | 20 | 7890
7 7 | 25 | 7890
8 8 | 34 | 7890

From Table 4 it is clear that $\Delta 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 5: Optimization of the TOLINTEG parameter (SHRINK = 2 2 and XLGRID). TOLINTEG |  $\Delta E$ ($\mathrm{K}$)
---|---
6 6 6 6 12 | 7889
7 7 7 7 14 | 7906
8 8 8 8 16 | 8063
9 9 9 9 18 | 7981
10 10 10 10 20 | 7890

As for the SHRINK parameter, the variations of TOLINTEG values does not
practically affect the final $\Delta E$ value. Consequently, its values have
been set equal to 7, 7, 7, 7, and 14.

Table 6: Optimization of the integration grid parameter (SHRINK = 2 2 and TOLINTEG = 7 7 7 7 14). GRID |  $\Delta E$ ($\mathrm{K}$)
---|---
XLGRID | 7951
XXLGRID | 7951

Again, the $\Delta E$ is practically unaffected by the changing the adopted
grid, 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 which interact with
each other just because of dispersive forces given the apolar nature of the
molecule. The results of Figure 15 clearly show that the B3LYP functional
needs some a posteriori correction (like the 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).

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
Angstrom.

###  A.5 ONIOM2 correction

Hereafter we summarize the single components (namely _BE_(CCSD(T), $\mu$-2H2O)
extrapolated, _BE_(DFT, $\mu$-2H2O), _BE_(DFT, ice) and finally the corrected
_BE_(ONIOM2)) of the ONIOM2 correction done with the extrapolation procedure
already described in § 2.

Table 7: Summary of the components used for the calculation of the ONIOM2 correction to the _BEs_. The last column shows the absolute value of the difference between the _BE_(DFT, $\mu$-2H2O) value and the corrected _BE_(ONIOM2). Species |  _BE_(CCSD(T), $\mu$-nH2O) extrapolated |  _BE_( DFT, $\mu$-nH2O) |  _BE_(DFT, ice) |  _BE_(ONIOM2) |  $|$Diff$|$
---|---|---|---|---|---
H2 | 843 | 796 | 1192 | 1240 | 47
CO | 1127 | 1165 | 2356 | 2318 | 38
CO2 | 1939 | 1817 | 3442 | 3564 | 121
OCS | 1354 | 1138 | 3478 | 3695 | 216
H2S | 2443 | 2875 | 5679 | 5248 | 432
HCN | 3300 | 3453 | 5795 | 5642 | 153
H2CO | 3199 | 3114 | 6490 | 6575 | 85
HCl | 3993 | 4475 | 6501 | 6019 | 482
NH3 | 4246 | 4389 | 7376 | 7233 | 143
CH3OH | 5196 | 5275 | 8681 | 8602 | 79
HCOOH | 6263 | 6223 | 9522 | 9562 | 40
HCONH2 | 5080 | 5056 | 9614 | 9638 | 24
OH• | 4492 | 4502 | 6542 | 6532 | 10
NH2• | 3021 | 2788 | 6043 | 6276 | 233
HCO• | 2339 | 2088 | 3474 | 3725 | 251

##  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 section. For some notable
cases, more than one initial structure has been modeled. All energetic
quantities are in Kelvin ($\mathrm{K}$), while distances are in Å. For every
molecule, the gas phase optimized geometry together with its ESP map are also
present. In the latter, the isovalue for the electron density is set equal to
$1\text{\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.

Fig. Set16. Adsorption on crystalline ice

”The complete figure set (24 images) is available in the online journal”.

Figure 16: I) Representation of the hydrogen (H2) molecule along with its ESP
surface, BSSE-corrected binding energies (_BE_) with (_disp_) and without (_no
disp_) dispersion. II) Top view of H2-(010) P-ice surface interaction (unit
cell highlighted in blue). III) Detail of the side view.

##  Appendix C Surface distortion energy on ASW

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

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 Species |  $\delta$Es crystalline ice |  $\langle$$\delta$ Es $\rangle$ amorphous ice
---|---|---
CH3OH | 662 | 2790
H2O | 565 | 3211
HCN | 589 | 1299
HCOOH | 878 | 1660
NH3 | 734 | 1082

##  Appendix D Basis set A-VTZ*

Table 9: Basis set for the H atom in CRYSTAL17 format used in this work H | 0 |  |
---|---|---|---
S | 3 | 1 |
|  | 34.06134 | 6.03E-03
5.123575 | 4.50E-02
1.164663 | 0.201897
S | 1 |  |
|  | 0.32723 | 1
S | 1 |  |
|  | 0.103072 | 1
P | 1 |  |
|  | 0.8 | 1
Table 10: Basis set for the C atom in CRYSTAL17 format used in this work C | 0 |  |
---|---|---|---
S | 5 | 1 |
|  | 8506.038 | 5.34E-04
1275.733 | 4.13E-03
290.3119 | 2.12E-02
82.0562 | 8.24E-02
26.47964 | 0.240129
S | 1 |  |
|  | 9.241459 | 1
S | 1 |  |
|  | 3.364353 | 1
S | 1 |  |
|  | 0.871742 | 1
S | 1 |  |
|  | 0.363524 | 1
S | 1 |  |
|  | 0.128731 | 1
P | 4 |  |
|  | 34.7095 | 5.33E-03
7.959088 | 3.59E-02
2.378697 | 0.142003
0.815401 | 0.342031
P | 1 |  |
|  | 0.289538 | 1
P | 1 |  |
|  | 0.100848 | 1
D | 1 |  |
|  | 1.6 | 1
D | 1 |  |
|  | 0.4 | 1
Table 11: Basis set for the N atom in CRYSTAL17 format used in this work N | 0 |  |
---|---|---|---
S | 5 | 1 |
|  | 11913.42 | -5.23E-04
1786.721 | -4.04E-03
406.5901 | -2.08E-02
114.9253 | -8.12E-02
37.10588 | -0.23871
S | 1 |  |
|  | 12.97168 | 1
S | 1 |  |
|  | 4.730229 | 1
S | 1 |  |
|  | 1.252518 | 1
S | 1 |  |
|  | 0.512601 | 1
S | 1 |  |
|  | 0.179397 | 1
P | 4 |  |
|  | 49.21876 | 5.55E-03
11.34894 | 3.81E-02
3.428509 | 0.149414
1.179951 | 0.348982
P | 1 |  |
|  | 0.417261 | 1
P | 1 |  |
|  | 0.142951 | 1
D | 1 |  |
|  | 2 | 1
D | 1 |  |
|  | 0.5 | 1
Table 12: Basis set for the O atom in CRYSTAL17 format used in this work O | 0 |  |
---|---|---|---
S | 5 | 1 |
|  | 15902.65 | 5.15E-04
2384.954 | 3.98E-03
542.7196 | 2.05E-02
153.4041 | 8.03E-02
49.54572 | 0.237668
S | 1 |  |
|  | 17.33965 | 1
S | 1 |  |
|  | 6.330336 | 1
S | 1 |  |
|  | 1.699588 | 1
S | 1 |  |
|  | 0.689545 | 1
S | 1 |  |
|  | 0.23936 | 1
P | 4 |  |
|  | 63.27052 | 6.07E-03
14.62331 | 4.19E-02
4.448952 | 0.161569
1.528151 | 0.356828
P | 1 |  |
|  | 0.529973 | 1
P | 1 |  |
|  | 0.175094 | 1
D | 1 |  |
|  | 2.4 | 1
D | 1 |  |
|  | 0.6 | 1
Table 13: Basis set for the S atom in CRYSTAL17 format used in this work S | 0 |  |
---|---|---|---
S | 5 |  |
|  | 103954 | 2.47E-04
15583.79 | 1.92E-03
3546.129 | 9.96E-03
1002.681 | 4.04E-02
324.9029 | 0.130675
S | 1 |  |
|  | 115.5123 | 1
S | 2 |  |
|  | 44.52821 | 0.504036
18.39789 | 0.230681
S | 1 |  |
|  | 5.510068 | 1
S | 1 |  |
|  | 2.125987 | 1
S | 1 |  |
|  | 0.436919 | 1
S | 1 |  |
|  | 0.157309 | 1
P | 5 |  |
|  | 606.6937 | 2.32E-03
143.507 | 1.86E-02
45.74616 | 8.60E-02
16.87291 | 0.252484
6.63992 | 0.446327
P | 1 |  |
|  | 2.672714 | 1
P | 1 |  |
|  | 1.000009 | 1
P | 1 |  |
|  | 0.354389 | 1
P | 1 |  |
|  | 0.116713 | 1
D | 1 |  |
|  | 1.1 | 1
D | 1 |  |
|  | 0.275 | 1
Table 14: Basis set for the Cl atom in CRYSTAL17 format used in this work Cl | 0 |  |
---|---|---|---
S | 5 |  |
|  | 117805.8 | 2.42E-04
17660.27 | 1.87E-03
4018.597 | 9.74E-03
1136.223 | 3.95E-02
368.1206 | 0.127972
S | 1 |  |
|  | 130.8615 | 1
S | 2 |  |
|  | 50.47901 | 0.428741
20.91681 | 0.196685
S | 1 |  |
|  | 6.353139 | 1
S | 1 |  |
|  | 2.494801 | 1
S | 1 |  |
|  | 0.543359 | 1
S | 1 |  |
|  | 0.194344 | 1
P | 5 |  |
|  | 681.0688 | 2.37E-03
161.1136 | 1.89E-02
51.38664 | 8.78E-02
18.95851 | 0.257074
3.003516 | 0.371524
P | 1 |  |
|  | 7.456529 | 1
P | 1 |  |
|  | 1.060936 | 1
P | 1 |  |
|  | 0.39452 | 1
P | 1 |  |
|  | 0.133233 | 1
D | 1 |  |
|  | 1.3 | 1
D | 1 |  |
|  | 0.325 | 1