{'fraction_non_alphanumeric': 0.03553497440896905, 'fraction_numerical': 0.015403363392639532, 'mean_word_length': 4.714763231197772, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Ab initio Calculation of Binding Energies of Interstellar Sulphur-Containing Species on Crystalline Water Ice Models

There are different environments in the interstellar medium (ISM), depending on the density, temperature and chemical composition. Among them, molecular clouds, often referred to as the cradle of stars, are paradigmatic environments relative to the chemical diversity and complexity in space. Indeed, there, radio to far-infrared observations revealed the presence of several molecules in the gas phase, while near-infrared spectroscopy detected the existence of submicron sized dust grains covered by H2O -dominated ice mantles. The interaction between gas-phase species and the surfaces of water ices is measured by the binding energy (BE), a crucial parameter in astrochemical modelling. In this work, the BEs of a set of sulphur-containing species on water ice mantles have been computed by adopting a periodic ab initio approach using a crystalline surface model. The Density Functional Theory (DFT)-based B3LYP-D3(BJ) functional was used for the prediction of the structures and energetics. DFT BEs were refined by adopting an ONIOM-like procedure to estimate them at CCSD(T) level toward complete basis set extrapolation, in which a very good correlation between values has been found. Moreover, we show that geometry optimization with the computationally cheaper HF-3c method followed by single point energy calculations at DFT to compute the BEs is a suitable cost-effective recipe to arrive at BE values of the same quality as those computed at full DFT level. Finally, computed data were compared with the available literature data.


Introduction
The interstellar medium (ISM) is the region between stars constituted by matter and radiation. The interstellar matter consists of either gaseous (atoms and molecules) or solid (dust grains of silicate or carbonaceous material) components. Several environments can be defined in the ISM according to their physical conditions, amongst them the dense molecular clouds, cold (10 K) and rarefied (10 3 cm −3 -10 7 cm −3 ) regions in which stars form. Because of the very low temperature, species heavier than He disappear from the gas phase due to freezeout (adsorption) on dust grain surfaces, forming a thick layer of ices. The ice composition is dominated by H 2 O ice, but also by small molecules like CO, CO 2 , CH 3 OH, CH 4 and NH 3 [4]. The significant presence of hydrogenated species on the ice mantles is due to the hydrogenation reactions taking place on the ice surfaces and involve atoms and small molecules that are adsorbed on grains [6].
Dust grains play a fundamental role in the chemistry of the ISM. Beyond their physical properties, like UV absorption and radio emission, they can help the occurrence of chemical reactions in three ways: i) they can concentrate chemical species on the surfaces, facilitating the encountering of the reactive species; ii) they can act as chemical catalysts, decreasing the activation energy of chemical reactions; and iii) they can be a third body, allowing the dissipation of the energy released by exothermic reactions without undermining the stability of the products.
Molecules on the grain surfaces require some energy to desorb. This energy can be thermal, chemical or radiative, but it has to exceed the interaction energy between the molecule and the surface since the evaporation rate is proportional to exp[−E B /(k T dust )], where T dust < 10 K and the binding energies are typically E B > 1000 K. If the amount of energy is a fraction of the binding energy, species can diffuse on the surfaces, which is one of the prerequisite for their subsequent reactivity [9]. Accordingly, one of the most important parameters in the description of the ISM are the binding energies (BEs), i.e., the interaction energy of chemical species on dust grains . Molecules adsorbed on grains can interact with them in two main ways: through rupture/formation of chemical bonds, so-called chemisorption, or through weak interactions, resulting in physisorption. This latter mechanism is the focus of this work and what our BEs are referred to. BE governs whether molecules remain stick on the grain surfaces as well as whether they diffuse and, therefore, react with other species. Thus, BEs hugely affect the gas and ices composition. Moreover, they are essential parameters to run astrochemical models at large aimed to understand the chemical composition of the ISM.
The importance of a careful estimation of BEs is thus evident. However, they have usually been estimated simply as a sum of the BEs of the single atoms or functional groups forming the molecule (in turn estimated from their polarizability when interacting with bare grain surfaces). In addition, the experimental procedure to collect BEs is not completely trustworthy, considering the impossibility to reproduce properly the ISM conditions in a laboratory chamber [9,18]. Ferrero et al. (2020) [11], have used quantum chemical computations as an alternative way to obtain accurate BE values for a set molecular species on water ices. The aim of this work is to enlarge the network of BEs, focusing on eight sulphur-containing molecules that play a major role in the dense clouds.

Methods
We adopted a periodic approach to model the adsorption of the sulphur-containing molecules (all of them being closed-shell species) on the ice surface. The hybrid B3LYP DFT functional [2,3,16], with the addition of Grimme's D3 empirical correction with Becke-Johnson (BJ) damping scheme to account for the dispersion interactions [13,14,19], in combination with the Ahlrichs triple zeta quality VTZ basis set, supplemented with a double set of polarization functions (Ahlrichs-VTZ*), was used for geometry optimizations and BE calculations (hereafter referred to as full B3LYP-D3(BJ) theory level).
Because of the high computational cost of DFT calculations, we also adopted the computationally cheaper semi-empirical HF-3c method in two ways [19]. This method is based on the Hartree-Fock (HF) energy computed with a minimal basis set, to which three empirical corrections (3c) are added. In the first case, HF-3c was used both in the optimization and the BE calculation (HF-3c//HF-3c theory level). In the second one, DFT single point energy calculations on the HF-3c optimized geometries were performed to compute BEs (B3LYP-D3(BJ)//HF-3c theory level). The latter approach was used to validate its performance in prevision of applying the same approach in amorphous ice models, which are prohibitively too expensive to be simulated at full DFT level.
Finally, to prove that DFT BEs are accurate enough, the single-and doubleelectronic excitation coupled-cluster method with an added perturbative description of triple excitations (CCSD(T)), in combination with a correlation consistent basis set extrapolated to complete basis set (CBS) limit, was used by applying a local correction at this level adopting an ONIOM2-correction approach [8], in which a small part of the system (i.e., the adsorbed molecule plus a few interacting water molecules) was treated at CCSD(T), while the rest at B3LYP-D3(BJ).
The BE is defined as the opposite of the interaction energy ∆E, which were corrected for the basis set superposition error (BSSE), due to using a finite basis set of localized gaussian functions.
All periodic calculations were performed with the crystal17 code [10], while CCSD(T) calculations were performed with the gaussian16 code [12].

Results
Interstellar ices are thought to be mostly formed by amorphous solid water [4]. We, however, have chosen a crystalline model for two main reasons: (i) crystalline structure are well defined due to symmetry constraints and computationally cheap, and (ii) no definite structure of an amorphous ice model is available. Nevertheless, regions rich in crystalline ices have been observed in protoplanetary disks and stellar outflows [17,20]. Our ice surface model derives from the bulk of the proton-ordered P-ice, which was cut along the (010) surface defining a 2D periodic slab model [5,22]. The thickness of the ice (10.9 Å) was chosen to converge its surface energy. The (010) P-ice supercell slab model consists of twelve atomic layers (24 water molecules) and its cell parameters are |a| = 8.980 Å and |b| = 7.081 Å. The structure of the ice is such as to ensure a null electric dipole along the non-periodic z-axis. This is a direct consequence of the symmetry of the system, which shows two identical faces both at the top and at the bottom of the model. Therefore, the adsorption was modeled only onto one face of the system.

BEs of the S-containing species
Here the BEs of 8 sulphur-containing species, i.e., H 2 S, CH 3 SH, H 2 S 2 , H 2 CS, C 3 S, SO 2 , CS and OCS, on the crystalline (010) water ice surface model were computed. Figures 1a and 1b show the optimized geometries of the computed complexes at full B3LYP-D3(BJ) level, while Table 1 reports the computed BE for these species.
According to the adsorption geometries, these species can be divided in two groups: i) molecules that can act as both hydrogen bond (H-bond) acceptors and donors, given the presence of an electronegative sulphur atom and one or more hydrogen atoms (H 2 S, CH 3 SH, H 2 S 2 and H 2 CS); and ii) molecules that can only act as H-bond acceptors, through an electronegative atom which is never the sulphur one (C 3 S, SO 2 , CS and OCS). We would like to highlight that this classification does not always correspond with the electrostatic (non dispersive) contribution to the BEs. Indeed, not only H-bond contributes to the definition of the electrostatic energy, but also the interaction between permanent dipoles and higher order multipoles. Additionally, dispersive forces are stronger as the volume of the molecule and the interaction surface with the ice increase, resulting in the articulated panorama appearing in Table 1. While the molecules both accepting and donating an H-bond show a BE spanning the limited range between 40 and 52 kJ mol −1 , the other group presents more variability. In fact, not only these molecules show very different dipoles, but they also differentiate in their mass and shape, giving rise to a dispersive force that can account for between 50% and 100% of the BE. The case of OCS is of particular interest because, by having the smallest dipole moment in the set of molecules, it is the less bound species to the ice. Not only it has the smallest binding energy, but also its interaction with the ice is almost completely (>90%) due to the dispersive forces.   Table 1: B3LYP-D3(BJ)/Ahlrichs-VTZ* computed BEs (in kJ mol −1 ) of the sulphur-containing species on the P-ice (010) (2x1) super cell model. The dispersive and non-dispersive contribution to the BE are also listed. Whenever more than one adsorption geometry was found, all BEs are listed and separated by a slash.

Molecule
Total

HF-3c method
As mentioned above, interstellar ice grains are mostly considered as amorphous solid water (ASW) systems. Unfortunately, ASW can only be simulated by requesting unit cells large enough to provide sufficient structural variability to mimic the amorphous nature. This increases the request of computer resources, making DFT methods almost impractical. Therefore, here we check our DFT results against those obtained with HF-3c, a much cheaper methodology to potentially be used when dealing with amorphous systems.
From the optimized DFT structures, we proceeded through the two HF-3c ways explained in the Methods section, i.e., at HF-3c//HF-3c and at B3LYP-D3(BJ)//HF-3c theory levels. The obtained results are summarized in Table 2. Figure 2a (left panel) shows the comparison between the BE computed at full B3LYP-D3(BJ) with the results computed at HF-3c//HF-3c. The correlation is only coarse, as it would be expected from the various approximations included in the HF-3c definition. In contrast, the performance of the B3LYP-D3(BJ)//HF-3c method is very good (see Figure 2a, right panel), showing its applicability to compute BEs on amorphous surface models in terms of computational cost and quality of the results, allowing moreover a wide exploration of different binding sites in an affordable way.

CCSD(T) correction
We finally performed some calculations at CCSD(T)/CBS (ONIOM2-correction, see Section 2) to check the accuracy of the B3LYP-D3(BJ) functional in describing the local interactions between the species and the ice and, at the same time, to refine the computed BEs. We chose to compute the correction only for some species, with the aim to span as much as possible the range of BEs and to sample geometries due to different contributions of dispersion forces. We are aware of the fact that five points are few for a meaningful statistics. However, when plotting the ONIOM2-corrected CCSD(T)/CBS BEs against those computed at the full B3LYP-D3(BJ) level, one can see an almost perfect correlation between the values, regardless of the nature of their interaction (see Figure 2b). Therefore, our results evidence the accuracy and reliability of the full B3LYP-D3(BJ) procedure adopted and, by extension, of the B3LYP-D3(BJ)//HF-3c scheme.

Discussion
Our computed BEs have been compared with the available literature data in Figure 3. In 2017, Penteado et al. [18] presented a list of recommended BEs, collecting data from previous works and adding an uncertainty range to each value. For H 2 S, OCS and SO 2 , the BEs were estimated from the experimental work of Collings et al. (2004) [7], while for all other species, data were based on the work of Hasegawa & Herbst (1993) [15], which collected BEs from previous works, including Allen & Robinson 1977 [1]. The uncertainty range for these values was set to half the BE value for species whose BE is less than 1000 K and to 500 K for all other cases. Comparison between these BE estimates with our computed BEs shows discrepancies, our values being systematically larger. This inconsistency should be taken with care, given the origin of the data published by Penteado et al. (2017). The cases of H 2 S, OCS and SO 2 , the ones extrapolated from thermal desorption process (TPD) experiments by Collings et al. (2004), are of interest. In the experiments, the species were deposited onto a pre-existing film of H 2 O at 8-10 K with an exposure that ensured a multilayer growth of each species. This is at variance with our assumption of a submonolayer adsorption. Therefore, while for OCS the difference amounts to about 30%, for H 2 S and SO 2 this percentage goes up to almost 60%. An additional reason explaining these large differences can be tracked down to the crystalline nature of our computer model compared to experimental data measured for AWS ice.  The paucity of different adsorption sites at the crystal ice and the enhanced hydrogen bond cooperativity compared with amorphous ices give BEs values higher than any experimental estimate. Additionally, the fact that our BE values do not account for zero-point energy (ZPE) corrections contributes to enlarge the discrepancies. The comparison with values based on the additive principle by Allen and Robinson in 1997 [1] present an average difference of 50%. Given the huge difference in the theoretical background standing behind these values, . [9,18,21] it is not straightforward explaining the differences. However, we can advance that the additive principle employed by Allen and Robinson [1] to estimate the BEs, although being a reasonable approach in 1977, it is currently not suitable to obtain BEs since the interaction of a molecule with a surface cannot be described as the sum of the interactions of the single atoms constituting the molecule.
In 2017, Wakelam et al. proposed a new chemical model to compute the binding energies [21]. Given the difficulty to reproduce a reliable structure of an amorphous ice, it was assumed that BEs on ASW are proportional to the interaction energy with one single water molecule, without ZPE nor BSSE correction. The approximate nature of this model was recovered by correcting the computed BEs through a scaling factor from the best fit with selected experimental BE values. The general trend by comparing our computed BE values with those of Wakelam et al. (2017) is a percentage difference of about 30%, the our values being larger. Although being corrected, it seems clear that using one water molecule to model an entire water ice surface is not a suitable choice since different important aspects such as H-bond cooperativity and surface long-range effect are neglected.
In 2018, Das et al. [9] simulated the water surface using four, and in some cases, six water molecules. For S-bearing species, only H 2 S and OCS were considered at MP2/aug-cc-pVDZ level on the water tetramer. Comparison of these values with our results with shows large differences. Reasons are twofold: i) the P-ice (010) surface exploits full hydrogen bond cooperativity, largely missing in the water tetramer; and ii) the adopted quantum mechanical methods are significantly different. Moreover BEs by Das et al. (2018) were not corrected for BSSE, which sure contribute to these discrepancies, nor for the ZPE [9].
The work of Ferrero et al. (2020) [11], presented for a set of chemical species a BE distribution by modeling their adsorption on both a crystalline and an amorphous surface ice models. Results showed that those on the crystalline ice were usually in the upper region of the range. Thus, we expect a similar behaviour for the adsorption of S-bearing species. Therefore, literature values should then be included in the range spanned by the distributions. With that being said, it is clear that the computation of BEs on an amorphous water ice model is necessary and already scheduled for the future, especially because these values have a pivotal role in the accuracy of the predictions provided by astrochemical models that deal with the chemical evolution in the ISM.

Conclusions
In this work we computed the BEs of 8 astrochemically-relevant S-containing species (i.e., H 2 S, CH 3 SH, H 2 S 2 , H 2 CS, C 3 S, SO 2 , CS and OCS) on water ice surfaces by means of quantum chemical calculations. We followed the same approach proposed by Ferrero et al. (2020) [11], in which the interstellar water ice was simulated through the (010) surface of the crystalline P-ice [5] as first step to subsequent studies in which amorphous ice systems will be employed. It is found that the S-containing species can be categorized in two groups: those acting as both H-bond acceptors and donors (H 2 S, CH 3 SH, H 2 S 2 , H 2 CS), and those acting only as H-bond acceptors (C 3 S, SO 2 , CS and OCS), in which the S atoms dos not participate in the interaction. Although a priori one could think that the former group would present larger BEs than the later ones, B3LYP-D3(BJ) results indicate that this is not the case since dispersive interactions are actually important, specially in the second group. CCSD(T) refinements on the BE values adopting an ONIOM2-correction approach indicate that those obtained at full B3LYP-D3(BJ) level are accurate enough. Similarly, comparison between B3LYP-D3(BJ)//HF-3c and full B3LYP-D3(BJ) BE values indicate that the former method is an optimum cost-effective method to compute BEs, in that case to be applicable in larger systems like the amorphous ones.
When comparing the computed BEs with experimental and other theoretical values from the literature we found major discrepancies. They arise from: i) the adoption of a crystalline ice model, which emphasizes the H-bond cooperativity compared with the amorphous ones, and ii) the different accounting for of long range order effects, highly present in crystalline systems while lost in the random network of the amorphous ice. Furthermore, the paucity of adsorption sites due to symmetry constraints as provided by the crystalline systems gives BE values that are at the upper extreme of the ranges compared with the larger surface morphological variability present in the amorphous ices. Such a constraint does not allow us to establish a statistical distribution of the BEs values, which is compulsory if a fair comparison with the experimental TPD data is foreseen. It is clear that the uncertainties on these values require more studies on the subject, especially because the binding energies have a pivotal role in the description of the chemistry occurring in the ISM.