  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 H\textsubscript{2}O-dominated
ices and dirtied by the presence of other volatile species. Whether a molecule
is in the gas or solid-phase depends on its binding energy (BE) on ice
surfaces. Thus, BEs are crucial parameters for the astrochemical models that
aim to reproduce the observed evolution of the ISM chemistry. In general, BEs
can be inferred either from experimental techniques or by theoretical
computations. In this work, we present a reliable computational methodology to
evaluate the BEs of a large set (21) of astrochemical relevant species. We
considered different periodic surface models of both crystalline and amorphous
nature to mimic the interstellar water ice mantles. Both models ensure that
hydrogen bond cooperativity is fully taken into account at variance with the
small ice cluster models. Density functional theory adopting both B3LYP-D3 and
M06-2X functionals was used to predict the species/ice structure and their BE.
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.
