Keywords

1 Introduction

The interstellar medium (ISM) is the region between the stars consisting of matter and energy. The ISM presents a vast range of physical conditions, with temperatures spanning from 10 K to 106 K, and densities varying from 10–4 to 108 cm–3 [1]. The interstellar matter is composed of gaseous material (such as atoms, molecules, and ions) and solid constituents (like submicron carbonaceous and silicate dust grains) [2]. In the cold (10 K) and dense (104 cm–3) environments of prestellar cores [3], also known as interstellar dark molecular clouds, volatile species such as H2O, CO, CO2, NH3, CH3OH, and CH4 freeze-out onto the dust grain surfaces forming ice mantles, predominantly of amorphous H2O ices and apolar CO ices [4].

As the gravitational collapse of the cold cloud proceeds, the interstellar gaseous material emits IR radiation, resulting in warm (100 K) inner envelopes around the low-mass protostar known as hot corinos, which provide viable conditions for rich grain surface chemistry, such as the formation of interstellar complex organic molecules (iCOMs), i.e., chemical compounds containing 6–13 atoms, of which at least one is C [5]. Although the complex species present in the cold regions are generally unsaturated and exotic, the warm cores in young stellar objects contain saturated (hydrogen-rich) complex organic molecules [6]. Despite being formed on the icy dust grains, as the temperature in the hot corinos reaches the mantle sublimation temperature, the species frozen in the ice mantles are released into the gas phase, where they can be detected and observed via their rotational spectra [7].

The molecules require energy to desorb from the grain surfaces. If the chemical, thermal, or radiative energy acquired by the molecule is higher than its interaction energy (i.e., binding energy, BE) with the icy surface, the species can evaporate into the gas phase. In case the acquired energy is a fraction of the BE, the molecules diffuse on the grain surface which is a critical step for chemical reactions to occur [7]. Therefore, obtaining highly accurate BE values of chemical compounds on dust grains is of utmost importance. Moreover, they are important input parameters in gas-grain astrochemical modelling for the characterization of the chemical composition and evolution of the ISM. The molecular species can interact with the icy grain surfaces by either forming chemical bonds (chemisorption), or via weaker (e.g., hydrogen bonding or dispersion) interactions (physisorption). This work focuses on the physical adsorption of iCOMs on a crystalline water ice surface model.

While techniques such as temperature programmed desorption (TPD) are used to estimate BEs experimentally, they are unable to accurately reproduce the extreme and diverse interstellar conditions in the laboratory setup. Hence, quantum chemical calculations provide a robust alternative to determine accurate BEs. This work, therefore, aims to add to the existing network of BEs those of some iCOMs recurrently found in astronomical observations, some of which, furthermore, might be regarded as potential prebiotic precursors.

2 Methods

In this work, the BEs of 13 closed-shell species have been calculated at various levels of theory adopting a periodic approach. Structural optimizations for the adsorption of iCOMs on the crystalline ice surface were carried out using the density functional theory (DFT), with the CRYSTAL17 program package [8]. The hybrid B3LYP functional [9,10,11] was utilized, along with Grimme’s DFT-D3 correction with the Becke-Johnson damping [12, 13] for the treatment of London-dispersion interactions, and an Ahlrichs triple zeta valence (TZV) basis set with polarization functions (henceforth referred to as B3LYP-D3(BJ)/Ahlrichs-TVZ*). Due to the use of finite basis set, the basis set superposition error (BSSE) was encountered, which has been treated using the counterpoise method (CP) [14] as follows:

$$ \Delta E^{CP} = \, \Delta E - BSSE $$
(1)

where ΔE is the interaction energy and ΔECP is the BSSE-corrected interaction energy using the CP method. The calculated binding energies are defined as the inverse of ΔECP:

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

In addition to the computationally costly B3LYP-D3(BJ)/Ahlrichs-TVZ* calculations, the significantly cheaper semi-empirical HF-3c treatment has also been applied for the geometry optimization of the adsorbate/ice complexes. This method uses the MINIX basis set, along with three a posteriori corrections pertaining to London dispersion interactions, BSSE, and short-range deficiencies arising due to the minimal basis set [15]. In order to further improve the HF-3c BEs, single-point energy calculations at B3LYP-D3(BJ)/Ahlrichs-TVZ* on the HF-3c optimized geometries were performed. This relatively accurate yet computationally cost-effective method (hereafter referred to as DFT//HF-3c) can be further employed to calculate BEs on more realistic (i.e., larger and amorphous) ice models, where full DFT treatment is unfeasibly costly.

In order to check the accuracy of the DFT calculated BEs, the CCSD(T) method has been employed, with the Gaussian16 software package [16], along with correlation consistent basis sets extrapolated to the complete basis set limit (CBS). Due to the exceedingly high computational cost of the CCSD(T)/CBS method, an ONIOM2-like approach has been applied, in which the molecules involved in the binding (the model zone, i.e., the iCOMs and three/four water molecules, depending on the size of the iCOM) are treated at CCSD(T), while computing at DFT for the whole system [17]. By following this ONIOM2-like scheme, the BEs are given as:

$$ \begin{gathered} BE\left( {CCSD\left( T \right), \, system} \right) = BE\left( {CCSD\left( T \right), \, model} \right) \hfill \\ + BE\left( {DFT, \, system} \right) - BE\left( {DFT, \, model} \right) \hfill \\ \end{gathered} $$
(3)

The water ice surface has been modelled using the bulk of the proton ordered crystalline structure of P-ice. A 2D-periodic slab model has been defined by cutting along the (010) surface [18,19,20]. The slab has a null electric dipole moment along the non-periodic z-axis. The slab model consists of twelve atomic layers involving 24 water molecules per unit cell, with cell parameters |a| = 8.980 Å and |b| = 7.081 Å. Due to the larger size of CH3CH2CN, a 2 × 1 supercell was used to ensure that the intermolecular interactions with the ice slab are predominantly considered, and spurious lateral interactions between adjacent unit cells are avoided.

3 Results

In this work, the BEs of 13 iCOMs (shown in Fig. 1) have been computed on the crystalline (010) water ice surface model. All the compounds are closed-shell species and, except HCOOH, comprise of 7–11 atoms. They contain different functional groups that determine their adsorption properties.

Fig. 1.
figure 1

Set of iCOMs used in this study for the calculation of their BEs on the crystalline water ice surface model.

Full size image

Based on likely hydrogen bond (H-bond) patterns, the iCOMs can be divided into two different categories: (i) iCOMs containing both H-bond donors and acceptors, due to the presence of hydrogen atoms chemically bonded to electronegative constituents, and (ii) iCOMs containing H-bond acceptors only, due to the presence of electronegative atoms such as O or N.

3.1 BEs of iCOMs at DFT

The B3LYP-D3(BJ)/Ahlrichs-TVZ* optimized geometries of the adsorbed species on the water ice surface are given in Fig. 2, while the calculated BEs and dispersion contribution are listed in Table 1. H-bonding and dispersion interactions are the two major binding driving forces. The general trend observed is that the species acting as both H-bond acceptors and donors have considerably higher BEs (around 80 kJ mol−1) than the molecules that can only act as H-bond acceptors (ranging 50–65 kJ mol−1). The formation of multiple H-bonds by different atoms of the adsorbate (namely, with one H-bond donor and at least one H-bond acceptor) leads to H-bond cooperativity, which strengthens the H-bond network in the system and, in turn, binds the iCOM more efficiently to the ice surface. The number of H-bonds significantly affects the values of the BEs, since for the complexes presenting three H-bonds, they are about 30 kJ mol−1 higher than those forming two H-bonds with the ice surface (e.g., CH3COOH and CH3CH2OH respectively). In this work, all the H-bonds between the adsorbate and the surface lie in the range of 1.65–2.16 Å.

Fig. 2.
figure 2figure 2

Optimized geometries of complex organic molecules adsorbed on crystalline (010) ice slab. The hydrogen bond lengths are given in angstrom. Color legend: Red, O; White, H; Grey, C; Blue, N (Color figures online). CH3CH2CN uses a 2 × 1 supercell due to the large system size.

Full size image
Table 1. Calculated BEs of iCOMs obtained at B3LYP-D3(BJ)/Ahlrichs-VTZ* level of theory. The contributions of dispersion and non-dispersion interactions are listed, along with the fraction of dispersion contribution to the total BE. The computed BEs are compared with BE values
Full size table

It is also worth noting that, in the species containing only H-bond acceptor groups, the contribution of dispersion interactions in the total BE is higher than those having acceptor and donor H-bond groups, which can form multiple H-bonds. Species such as HCOOH and HOCH2CHO have significantly lower extent of dispersion (37% and 36% of the total respective BEs), whereas dispersion plays the main role (more than the non-dispersive interactions) in molecules such as CH3CH2CHO and CH3COOCH3 (53% and 57% of the total respective BEs). The species containing only H-bond acceptor groups exhibit a more non-polar character than the species containing both H-bond donor and acceptor groups, because the former are CH-containing compounds with heteroatoms. Accordingly, the dispersive contribution to the binding energy of this set of compounds is higher than the H-bond donor and acceptor groups, while for the species in the latter group, the main contributor to the binding energy is H-bonding.

Our computed BE values have been compared with data available in the literature. There are some discrepancies between the previously reported and our computed values. They can mainly be attributed to the difference of the water ice surface models employed. In this work, we have employed a periodic crystalline ice slab consisting of 24 H2O molecules per unit cell. A similar periodic strategy was used in Lattelais et al. [21]. However, they used a different quantum mechanical methodology: calculations at GGA PW91 level of theory using planewaves as basis set. The main difference between the two methods is that Lattelais’ calculations lack the inclusion of dispersive forces, which have been seen in this work to be crucial. Consequently, in most of the cases, their values are lower than ours. In Wakelam et al. [22], M06-2X was employed, which is a meta-hybrid DFT method and includes dispersive forces in their definition. However, the main difference is in the water ice surface model: they used a single water molecule. Accordingly, long-range effects such as H-bond cooperativity and extended, actual iCOM-surface dispersion interactions are neglected. Thus, Wakelam’s values are also lower than those computed by us.

3.2 Comparison with HF-3c Method

Since most of the water ice in the ISM exists in the amorphous form [4], it is imperative to compute the BEs of the species of interest on the amorphous water ice system. Due to the lack of definite long-range order, an amorphous water ice model would require a much larger unit cell to accurately simulate the structural variability of the amorphous ice mantle. In such cases, methods like DFT would be unfeasible and, accordingly, comparatively cheaper methods such as HF-3c are more practical. Thus, we test the accuracy of the BEs provided by the DFT//HF-3c scheme against those from B3LYP-D3(BJ)/Ahlrichs-TVZ* (here referred to as DFT//DFT). Computed BEs are listed in Table 2.

Table 2. Computed BEs of iCOMs at different levels of theory. The energies are given in kJ mol−1. The values marked with asterisks present BEs on a different HF-3c optimized structure (more details in the text).
Full size table

Interestingly, in some cases, HF-3c optimized structures differed from the DFT optimized ones, hence resulting in inconsistencies in the computed BEs with both methods. For these cases, we used the DFT optimized structures as the initial guess structure for a new geometry optimization at HF-3c, thus obtaining another BE value that correlates better with the DFT one. Table 2 presents both BEs, with the different geometries marked with asterisks.

3.3 CCSD(T) Refinement

Fig. 3.
figure 3

Correlation among DFT//HF-3c, DFT//DFT, and CCSD(T)//DFT BEs. The BEs are given in kJ mol−1.

Full size image

Table 2 also presents the BEs computed with the ONIOM2-like CCSD(T)/CBS scheme (referred to as CCSD(T)//DFT) to check the accuracy of the DFT//DFT BEs. Our ONIOM2-like approach includes the effects of the whole system at B3LYP-D3 level of theory and provides energies comparable to a full CCSD(T) treatment. Figure 3 shows the correlation among the computed BE values with the different methodologies: DFT//DFT, DFT//HF-3c, and ONIOM2-like CCSD(T)//DFT. Results indicate that the CCSD(T)/CBS computed BEs are in an excellent correlation with the DFT//DFT values, and the cost-effective DFT//HF-3c BEs are in good correlation with DFT//DFT and CCSD(T)//DFT too. Hence, it turns out that DFT//HF-3c is a reliable method for the calculation of BEs of iCOMs on large amorphous water ice surfaces.

4 Conclusion

In this work, the BEs of 13 iCOMs (i.e., CH3COOH, CH3CH2OH, HOCH2CH2OH, HCOOH, HOCH2CHO, c-CH2CH2O, CH3CHO, CH3COCH3, CH3CH2CN, CH3COOCH3, HCOOCH3, CH3NCO, and CH3CH2CHO) on a P-ice based crystalline (010) water ice surface model have been computed by means of periodic quantum chemical calculations. With the aim to provide accurate BEs to be used as input parameters in astrochemical models, the B3LYP-D3(BJ)/Ahlrichs-TVZ* level of theory has been employed. Additionally, a computationally cheaper DFT//HF-3c scheme has been used in order to compare its accuracy with the DFT level for further application in amorphous water ice systems. BE refinements have been carried out at CCSD(T)/CBS level of theory using the ONIOM2 methodology to compare the DFT results with the “gold-standard” accuracy provided by the coupled cluster method.

The results reveal that some iCOMs (CH3COOH, CH3CH2OH, HOCH2CH2OH, HCOOH, and HOCH2CHO) can act as both H-bond donors and acceptors when adsorbed on the water ice surface. The BEs of such species are significantly higher (nearly 25–30 kJ mol−1) than the molecules that can only act as H-bond acceptors. This indicates that the former species are much strongly bonded to the ice surface than the latter counterparts due to larger H-bond cooperativity. On the other hand, it has been observed that dispersion accounts for a major fraction of the total BE in the case of iCOMs with H-bond acceptors only, while the dispersion contribution is lower in species with high BE values.

The comparison of BEs between those obtained at the full DFT level with those refined at CCSDT//DFT shows an excellent correlation, indicating that the B3LYP-D3(BJ) functional yields accurate energy estimation for the adsorbate/ice systems. Furthermore, DFT//HF-3c results correlate well with the full DFT treated ones, thereby indicating that the DFT//HF-3c scheme is a cost-effective level of theory to be utilized for the calculation of BEs on amorphous ice systems, which require a larger unit cell to account for the lack of long-range order and varied binding sites on the surface.