  The knowledge of the binding energy of molecules on astrophysically relevant
ices can help to obtain an estimate of the desorption rate, i.e. the molecules
residence time on the surface. This represents an important parameter for
astrochemical models, crucial to determine the chemical fate of complex organic
molecules formed on dust grains and observed in the densest regions of the
interstellar medium. In this work, we propose a new robust procedure to study
the interaction of atoms and molecules with interstellar ices, based on
\textit{ab initio} molecular dynamics and density functional theory, validated
by high-level \textit{ab initio} methods at a CCSD(T)/CBS level. We have
applied this procedure to hydrogen fluoride (HF), a promising tracer of the
molecular content of galaxies. In total we found 13 unique equilibrium
structures of HF binding to small water clusters of up to 4 molecules, with
binding energies ranging from 1208 to 7162 K. We computed a 22-molecules model
of amorphous solid water (ASW) surface using \textit{ab initio} molecular
dynamics simulations and carried out a systematic analysis of the binding sites
of HF, in terms of binding modes and binding energies. Considering 10 different
water clusters, we found a binding energy distribution with an average value of
$5313\pm74$ K, and a dispersion of $921\pm115$ K. Finally, the effect of the
electrostatic field of the 22 water molecules on the binding energies was
investigated incrementally by symmetry adapted perturbation theory, in order to
gauge the effect of the water environment. The results indicate that the extent
of the electrostatic interaction of HF with ASW depends strongly on the
properties of the binding site. We expect that this work will provide a solid
foundation for a systematic development of a binding energy distribution
database of molecules on interstellar surfaces.
