key: cord-0058689-akewuhl1
authors: Martínez-Bachs, Berta; Ferrero, Stefano; Rimola, Albert
title: Binding Energies of N-Bearing Astrochemically-Relevant Molecules on Water Interstellar Ice Models. A Computational Study
date: 2020-08-20
journal: Computational Science and Its Applications - ICCSA 2020
DOI: 10.1007/978-3-030-58808-3_49
sha: 458cd14379a72b8eac4c224dd3a69fff06179619
doc_id: 58689
cord_uid: akewuhl1

Different molecules and radical species have been detected and identified in the interstellar medium (ISM) through rotational and infrared techniques. The presence of dust grains in ISM have also been observed, which in the denser regions are covered by thick ice mantles consisting mainly of water molecules. Quantifying the interaction of gas phase species with water ice is essential to accurately model the evolution of the ISM. Hence, the importance to obtain accurate binding energies. In this contribution, the binding energies (BEs) of different N-containing interstellar species on water ice surfaces are presented. They have been calculated by means of reliable periodic quantum chemical simulations. Two different methods have been used: i) hybrid DFT methods, i.e., B3LYP for closed-shells species and M06-2X for open-shell radicals, and ii) the cost-effective DFT//HF-3C, i.e., geometry optimization at HF-3c followed by a single point energy refinement at the DFT level. As a first step, BEs have been computed on a crystalline water ice surface model of a thickness of 12 atomic layers, in which different super cells have been adopted (involving 36 or 72 water molecules per unit cell), depending on the size of the adsorbate. Moreover, an ONIOM-like local approximation combining CCSD(T) with DFT functional as high and low theory levels, respectively, has been adopted to correct BEs.

Interstellar medium (ISM) is the matter and radiation present in the regions between stars. Interstellar matter can be found as ionic, atomic and molecular gas forms but also as tiny dust grains. Interstellar matter is not homogeneous and its temperature ranges between 10 K-10 6 K and atomic densities between 10 −4 cm −3 -10 8 cm −3 depending on the environment and region [1] .

In the coldest and densest regions of the ISM, cloud-like environments are found. These are regions of the Universe with extension between 100 parsecs and 300 light years, in which gaseous atoms and molecules and dust grains are mixed together [2] . In these regions, dust grains are covered by thick ice mantles, which are mainly composed by water (H 2 O) but also by other minor volatile species, such as carbon monoxide (CO), carbon dioxide (CO 2 ), ammonia (NH 3 ), methane (CH 4 ) and methanol (CH 3 OH) [3] .

Infrared-based astronomical observations reveal that these ices are in an amorphous structural state. They consist of different ice mantles coating the core of the dust grains. Ice mantles can be in the form of (i) water-rich, so-called polar phases, which mainly contain H 2 O but also CH 3 OH and NH 3 and where hydrogenation reactions dominate, or (ii) water-poor mantles, so-called apolar phases, which contain mainly CO but also CO 2 and HCOOH and where the reactions involving CO are dominating. Furthermore, molecules from the gas phase can adsorb on the ice surfaces, diffuse on them and react to form other, more complex species [3] .

Astrochemical models are computational tools the purpose of which is to reproduce and represent the physical and chemical behaviour of interstellar clouds. Initially, astrochemical models were built to represent the abundances of early gas phase species [4, 5] . However, after noticing the fundamental role of dust grains in the increase of the molecular complexity in the ISM, novel astrochemical models have been developed that account for the chemistry on the dust grains [6] [7] [8] .

To obtain astrochemical modelling results, a set of physical and energetic parameters must be introduced as input data. Binding energy (BE) values are input parameters of paramount importance in the gas-grain astrochemical models as they dictate ice accretion, surface diffusion and desorption events, which are fundamental elementary physico-chemical processes in the grain surface chemistry. One of the major current problems is the lack of accurate molecule/ice BE values [8] . Temperature programmed desorption is the usual technique to determine BEs experimentally, but laboratory experiments present severe limitations as they are not capable to reproduce neither the actual interstellar conditions nor realistic ice grain analogues [10] . This leads to obtain unreliable BEs, which is a dramatic drawback in astrochemical modelling, and hence the importance to accurately quantify the interaction between gas-phase species and ice surfaces.

In this work, results on the binding energies (BEs) of different N-containing species of interstellar interest on water ice surfaces modelled by a crystalline system adopting a periodic slab approach are presented.

BEs have been calculated by means of robust periodic quantum chemical simulations adopting different methods and density functional theory (DFT) functionals. One quantum chemical method is the novel cost-effective HF-3c, in which the electronic energy is grounded on the Hartree-Fock framework, but three corrections are added in a posteriori way [11] . Two DFT functionals have also been used: B3LYP for closed-shell species and M06-2X for radicals open-shell species. Both methods are hybrid DFT functionals as they incorporate a certain amount of exact exchange, 20% and 45%, respectively [12] [13] [14] . For B3LYP, the D3 empirical correction has been added to account for dispersive forces [15] . For both DFT functionals, the Ahlrichs TZV basis with polarization functions were used, while for HF-3c the MINIX basis set was used. HF-3c values do not suffer from the basis set superposition error (BSSE) but DFT values do. Accordingly, for these cases, single point energy calculations to obtain BSSE-corrected BE values have been performed [16] .

Truncated Coupled Cluster (CC) methods are one of the most accurate electronic structure methods but due to its high computational cost they can only be applied for small chemical systems. In this work the CCSD(T), a method that treats the single and double excitations with a CC formalism and the triple ones perturbatively, has also been used to refine the DFT BEs adopting an ONIOM-like approach [17] , hereafter referred to as O[CCSD(T)//DFT]. In this apporach, a molecular system including the adsorbate species and the binding site region has been extracted from the periodic system and calculated at CCSD(T). The O[CCSD(T)//DFT] BEs are obtained as:

where the model system is the molecular "adsorbate+binding site" region and the real system is the actual periodic system, and the high method is CCSD(T) and the low method the DFT ones. For these calculations, we used the DFT optimized geometries. The surface mimicking the icy water dust grain has been modelled adopting a periodic approach. The water ice surface is based on the crystalline P-ice structure. P-ice is the proton ordered analogue to the hexagonal (proton-disordered) water ice, which has been demonstrated to reproduce fairly well the physico-chemical features of crystalline water ice [18] . The surface model has been generated by cutting out the P-ice 3D periodic bulk structure perpendicular to the [010] direction, resulting with the (010) slab surface model. We chose the (010) slab model because it is one of the most stable planes of the P-ice and because it has no dipole component perpendicular to them [19] . The slab model consists of twelve atomic layers and consists of 72 atoms (24 water molecules) per unit cell. It is worth mentioning that a structurally proper model based on the hexagonal ice (Ih) requires a large unit cell as all the proton arrangements must be covered. Indeed, as mentioned above Ih is a proton-disordered system. Since the position of the protons in the ice follows the Bernal-Fowler-Pauling (BFP) rules [20, 21] , i.e., i) each water molecule is oriented such that its two hydrogen atoms are directed approximately toward two of the four surrounding oxygen atoms arranged almost in a tetrahedron, ii) only one hydrogen atom is present on each OÁÁÁO linkage and iii) each oxygen atom has two nearest neighbouring hydrogen atoms such that the water molecule structure is preserved, a structural model for Ih envisaging all the hypothetical proton-ordered configurations accounting the BFP rules would be very large and, accordingly, costly to treat at a quantum mechanical level. Because of that, using P-ice is a good compromise between accuracy and computational cost.

The BEs of eight N-containing species, i.e., HNO, N 2 O, NO, NO 2 , N, NH, NH 2 and HNC, on the water crystalline slab model have been computed. These species can be classified in two groups, according to their capability to establish hydrogen bonds (Hbonds): molecules that can only act as H-bond acceptors (N, NO, NO 2 and N 2 O, hereafter referred to as Hb-A), and molecules that can act as both H-bond donors and acceptors (HNO, NH, NH 2 and HNC, hereafter referred to as Hb-AD).

As mentioned above, three different quantum chemical methodologies have been used: i) a full DFT treatment, using B3LYP for closed-shell species and M06-2X functional for open-shell species, ii) the cost-effective DFT//HF-3c treatment, in which geometry optimizations were carried out at a HF-3c level followed by a single point energy calculation at the corresponding DFT level, and iii) an ONIOM-like local approximation with CCSD(T) as high level of theory and the DFT functionals as low level of theory, from here on referred as O[CCSD(T)//DFT]. The purpose is to obtain accurate and robust BE values as well as to carry out a calibration study to assess the accuracy of the cost-effective DFT//HF-3c method, that is, if the computed BE values are in fair agreement with the full DFT (DFT//DFT) ones and the O[CCSD(T)//DFT] ones.

The structures of the complexes formed by adsorption of the Hb-A and Hb-AD species on the water ice surface model are shown in Figs. 1 and 2 , respectively. Table 1 reports the BE values obtained at both DFT//DFT, DFT//HF-3c and O[CCSD(T)//DFT] methods.

As expected, for all cases, H-bond interactions and the dispersion forces are the main intermolecular forces that dictate the adsorption of the species on the water ice surface. Optimized complexes at full DFT theory and at the HF-3c method present differences in the H-bond distance, which in most of the cases lay within 0.2 Å . In addition to the intrinsic differences in the definition of each method, the fact that HF-3c does not suffer from BSSE in the optimization processes while this is indeed the case for the DFT methods also has an influence in the observed H-bond differences.

In relation to the calculated BE values, the general trend is that Hb-AD species present larger BEs than Hb-A species. Taking the O[CCSD(T)//DFT] BE values, the former group have BEs spanning the 70.4-35.1 kJ mol −1 range while the latter group the 29.2-13.8 kJ mol −1 range. This is because the Hb-AD species can establish more and more efficient H-bond interactions with the surface than the Hb-A ones. More Hbonds because for all the Hb-AD cases two H-bonds are established, i.e., one as Hbond donor and one as H-bond acceptor. This is not the case for Hb-A species, as they can only establish one or two H-bonds, in this later case the acceptor atom being the same. And more efficient H-bonds because the fact to act as both H-bond donor and acceptor enables H-bond cooperativity, which reinforces the H-bond network in which the adsorbate species is involved. This is reflected by the H-bond distances identified in both groups of species: Hb-AD species establish H-bond interactions of between 1.6-2.0 Å, while Hb-A species of between 2.2-3.0 Å (at HF-3c theory level), meaning that this later group present weak H-bonds and accordingly dispersion is a major contributor to the BEs. In the Hb-AD group, HNC is the species presenting the largest BE value because it presents the strongest H-bond interaction (1.59 Å at HF-3c). In contrast, NH is the species with the lowest Hb-AD BE because their H-bonds are relatively weak (about 2.0 Å at HF-3c). In the Hb-A group, NO 2 presents the largest BE value because, being one of the largest molecules (and hence presenting large dispersion contributions), it exhibits the "strongest" H-bond (2.24 Å) among the species of its group. In contrast, N has the lowest BE because the H-bonds are actually very weak (2.77 and 2.98 Å at HF-3c) and dispersion contribution is very minor as it is a single atom.

An identified general trend is that DFT//DFT BE values are slightly more favourable than the DFT/HF-3c ones. This is due to the effects of the BSSE in the optimization at either full DFT or HF-3c theory levels. As mentioned above, optimized HF-3c geometries, which do not suffer from BSSE, generally present H-bond distances larger than the DFT optimized ones. Because of these larger H-bonds in the HF-3c geometries, the DFT single point energies provide lower BE energies than the DFT// DFT ones. Despite this fact, BE values obtained with the two methods are in almost perfect correlation, as one can see in Fig. 3 . Moreover, both methods also correlate very well with the BE values obtained through the O[CCSD(T)//DFT] approach (see Fig. 3 ), allowing us to conclude that the cost-effective DFT//HF-3c is a suitable method to obtain robust and accurate enough BE values. As a final interesting comment, we would like to highlight that B3LYP and M06-2X functional methods, based on the comparison with the O[CCSD(T)//DFT] approach, describe well the electronic structure of respective closed-shell and open-shell systems. 

In this contribution, results on the binding energy (BE) values of eight N-containing species of astrochemical relevance (i.e., N, NO, NO 2 , N 2 O, NH, HNO, NH 2 , HNC) on water icy grain surfaces, in this case modelled by a periodic (010) crystalline water ice slab model of P-ice, are presented. The main purpose is to assess the robustness of the cost-effective DFT//HF-3c approach (i.e., optimization at the composite HF-3c theory level followed by a single point energy calculation at DFT theory level) to calculate accurate BEs. To this end, DFT//HF-3c BEs have been compared to those provided at DFT//DFT as well as to those obtained through a CCSD(T) correction adopting an O [CCSD(T)//DFT] approach. The N-containing molecules considered here, in their adsorption on the water ice surface can act either as just H-bond acceptors (Hb-A, N, NO, NO 2 and N 2 O) or as Hbond acceptors and donors (Hb-AD, NH, HNO, NH 2 and HNC). Consequently, the latter group presents more favorable BE values than the former one. Indeed, while Hb-AD species establish relatively strong H-bond interactions with the surfaces, in Hb-Abased complexes H-bond interactions are relatively weak, leaving the major contribution to the total BE in the dispersive forces.

Despite the small differences between the optimized geometries at HF-3c and at full DFT levels, BE energies computed at DFT//HF-3c and at DFT//DFT are in almost perfect agreement. Moreover, BE values obtained with the two methods also agree very well with those obtained adopting the O[CCSD(T)//DFT] approach. These results indicate that DFT//HF-3c is a reliable cost-effective method to obtain accurate BEs, and accordingly it can be used in more realistic interstellar water ices, such as amorphous surfaces, which require the use of larger unit cells since the complex surface heterogeneity and binding sites diversities have to be accounted for.

Introduction: astrochemistry

Complex organic interstellar molecules

Astrochemistry of dust, ice and gas: introduction and overview

The density of molecules in interstellar space

The formation and depletion of molecules in dense interstellar clouds

Molecular dynamics simulation of the H2 recombination on a graphite surface

Models of gas-grain chemistry in dense interstellar clouds with complex organic molecules

A three-phase chemical model of hot cores: the formation of glycine

Grain surface models and data for astrochemistry

Corrected small basis set Hartree-Fock method for large systems

Density-functional exchange-energy approximation with correct asymptotic behavior

A new mixing of Hartree-Fock and local density-functional theories

The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other function

Accurate description of van der Waals complexes by density functional theory including empirical corrections

The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors

Density Functional Theory in Quantum Chemistry

Proton-ordered ice structures at zero pressure. A quantum-mechanical investigation

IR spectral fingerprint of carbon monoxide in interstellar water-ice models

A theory of water and ionic solution, with particular reference to hydrogen and hydroxyl ions

The structure and entropy of ice and of other crystals with some randomness of atomic arrangement