Gliding on Ice in Search of Accurate and Cost-Effective Computational Methods for Astrochemistry on Grains: The Puzzling Case of the HCN Isomerization.

The isomerization of hydrogen cyanide to hydrogen isocyanide on icy grain surfaces is investigated by an accurate composite method (jun-Cheap) rooted in the coupled cluster ansatz and by density functional approaches. After benchmarking density functional predictions of both geometries and reaction energies against jun-Cheap results for the relatively small model system HCN···(H2O)2, the best performing DFT methods are selected. A large cluster containing 20 water molecules is then employed within a QM/QM' approach to include a realistic environment mimicking the surface of icy grains. Our results indicate that four water molecules are directly involved in a proton relay mechanism, which strongly reduces the activation energy with respect to the direct hydrogen transfer occurring in the isolated molecule. Further extension of the size of the cluster up to 192 water molecules in the framework of a three-layer QM/QM'/MM model has a negligible effect on the energy barrier ruling the isomerization. Computation of reaction rates by the transition state theory indicates that on icy surfaces, the isomerization of HNC to HCN could occur quite easily even at low temperatures thanks to the reduced activation energy that can be effectively overcome by tunneling.

Introduction

In 2018, McGuire published a census of Interstellar, Circumstellar, Extragalactic, Protoplanetary Disks, and Exoplanetary Molecules[1] including more than 200 molecules (containing from 2 to 70 atoms) and this number is steadily increasing thanks to the modern technologies of new observatory telescopes.[2] The identification of many interstellar complex organic molecules (iCOMs) defeated the old and general idea that the interstellar medium (ISM) was an empty vial where chemical reactivity could not operate. Questions about the formation of iCOMs in such extreme conditions and the evolution of molecular complexity fueled the curiosity of astrochemists all over the world.[3] While gas-phase reactions seemed the obvious choice to explore the formation pathways of molecular systems in such rarefied environments, the ubiquitous presence of dust and grains and the mismatch between some observations and the molecular abundances predicted by gas-phase models have boosted the role of solid-state chemistry.[4,5] Since the discovery of the catalytic role of grains for H2 formation,[6,7] astrochemists and physicists have struggled looking for gas-grain models that could provide a comprehensive picture of chemical processes in the ISM. At the low temperatures of molecular clouds, molecules in the gas phase accrete icy mantles freezing out onto grain surfaces[8,9] and leading to porous and amorphous icy surfaces,[10−12] which can host local reactants triggering a molecular reactivity not feasible in the gas phase. The composition and morphological features make the simulation of these icy structures a great challenge in this field.[9,13]

The difficulty of performing experimental studies for systems capable of mimicking the harsh conditions of the ISM, calls for computational simulations of periodic surfaces and/or suitable model clusters able to take into the proper account the main structural features responsible for the chemistry at the interface.[14,15] This translates into the necessity of simulating extended systems, thus making the computational burden prohibitive for the accurate state-of-the-art methods developed for isolated molecules.[16] Because water is the main component of polar icy mantles,[17,18] a lot of efforts have been devoted to the investigation of the adsorption and formation of iCOMs on water clusters used to mimic interstellar ices. The structures of H2O clusters containing up to 22 atoms have been worked out from molecular dynamics simulations and made available in online databases.[19] Some years ago, Rimola et al. studied iCOM formation pathways on clusters including up to 33 water molecules obtained by combining two (H2O)18 clusters taken from the (010) surface of ice-XI[20] and removing three molecules to facilitate the construction of the final cluster.[21] Furthermore, attempts to include the structural modifications induced by UV and cosmic rays photo-processing have been made by means of small radical and ionized water clusters.[22] More recently, molecular dynamics has been used to model amorphous water ices[23] and to simulate mixed CO/H2O ices.[24] Adsorption energies on clusters of larger size have been evaluated by a two-layer our own N-layered integrated molecular orbital molecular mechanics (ONIOM) model, with the higher-level layer treated by means of density functional theory (DFT), and the lower-level one described through molecular mechanics (MM)[25,26] or semiempirical quantum chemical methods.[27]

While coupled cluster theory including full treatment of single and double excitations together with perturbative estimation of triple excitations [CCSD(T)], possibly in conjunction with composite schemes to estimate the complete basis set (CBS) limit, is considered the gold-standard for accurate predictions,[28] the size of the systems to be dealt with in the case of ice-mediated chemistry makes density functional theory the only viable route in terms of accuracy to computational cost trade-off. As is well known, the reliability of DFT strongly depends on the specific system and properties at hand and on the choice of the density functional (DF) among an ever increasing number of possible formulations. In this respect, benchmark is a fundamental step for ranking the reliability of DFT model chemistries, also in connection with the computational cost, and hence it represents a very active field of research.

Concerning the specific topic of adsorption and reactivity of iCOMs on interstellar ice analogues, to the best of our knowledge, systematic benchmark studies are still lacking. In this connection, Enrique-Romero et al.[29] performed a calibration analysis of radical–water interactions and activation energy for NH2 + HCO and CH3 + HCO reactions in the presence of one and two water molecules. They tested the accuracy of B3LYP and BHLYP functionals (both with and without dispersion corrections) in conjunction with the 6-311++G(2df,2pd) basis set taking CASPT2/cc-pVTZ and CCSD(T)/aug-cc-pVTZ levels of theory as reference. That analysis was focused on the interaction and activation energies, while recent works have highlighted that reliable geometries are fundamental prerequisites for accurate thermochemistry and kinetics.[16] In this respect, the B3LYP functional can be unable to predict correct structures for van der Waals complexes[30] and transition states.[31] Furthermore, the use of CCSD(T)/triple-ζ energies cannot be recommended as a reference in benchmark studies because basis set truncation and lack of core–valence correlation limit the accuracy, thus introducing a bias in the reference values. This issue can be overcome by resorting to composite methods that aim at minimizing the errors relying on well-tested additive approximations.[16,32]

In this work, we assess the performances of several DFT model chemistries in evaluating the structural and energetic aspects of ice-mediated interstellar reactions employing the HCN ⇌ HNC isomerization catalyzed by water molecules as a paradigmatic process. On the one side, this can be considered a model for more complex reactions mediated by ice surfaces and, on the other side, the chosen system is small enough to allow the exploitation of state-of-the-art composite methods to generate accurate reference values for both geometries and reaction energies. The HCN ⇌ HNC isomerization has been widely studied because the observed HNC/HCN ratio in the ISM cannot be predicted on the basis of the proposed gas-phase mechanisms. Moreover, both HCN and HNC can be involved in the formation of amino acid precursors in the Strecker synthesis of glycine.[33,34] Gardebien and Sevin investigated the process for the isolated molecule and with explicit inclusion of two to four water molecules[35] finding that the most favorable mechanism consists of a one-step path involving a proton relay mediated by the water cluster. Koch et al. employed a more realistic model including seven additional water molecules to simulate the local environment of the icy surface and employing the polarizable continuum model (PCM) to account for bulk effects.[36] According to the available data, the water cluster acts as a catalyst lowering the energy barrier with respect to the gas phase, an effect that progressively smooths increasing the number of H2O molecules. Intermolecular hydrogen bond drives both the interaction of HCN and HNC with the ice surface and the isomerization process. This represents the most common mechanism through which molecules adsorb and react on ISM polar ices.

On these grounds, we decided to perform a detailed study of the HCN ⇌ HNC isomerization by state-of-the-art quantum chemical methods and realistic cluster models. The work is organized as follows: the computational methods are described in Section , while the outcomes of the benchmark are detailed in Section concerning both geometries and energies, thus leading to the identification of the best performing DFT model chemistries in terms of the trade-off between accuracy and computational cost. Despite the fact that the benchmark is carried out on a simplified model, the outcomes are expected to be of general validity, especially with respect to the relative performances of the tested methods which can then be transferred to larger H2O clusters. With this in mind, at the end of Section , the best performing methods are employed to simulate the HCN ⇌ HNC isomerization catalyzed by a cluster of 20 water molecules and then further embedded in a 172 water slab described through MM. Finally, reaction rates are computed in the framework of the transition state theory (TST) including tunneling.

Computational Methodology

For the benchmark study, we selected 10 DFs belonging to different families: two hybrids (B3LYP, BHLYP),[37−39] a long-range corrected DF (ωB97X-D),[40] three meta-hybrids (PW6B95,[41] BMK,[42] and M06-2X[43]), one meta-NGA (MN15[44]), the B2PLYP,[45] and the two spin-component-scaled (DSD-PBEP86 and revDSD-PBEP86)[46,47] double hybrids. To test the accuracy to computational cost trade-off, for each functional, six basis sets have been considered. In particular, we selected Dunning’s aug-cc-pVnZ basis sets (n = D, T)[48,49] as well as the corresponding jun- and jul-modifications from Truhlar’s calendar family.[50] All the DFT calculations include empirical dispersion corrections according to the DFT-D3 scheme proposed by Grimme[51] with the Becke-Johnson damping function (BJ),[52,53] which are fundamental for the correct prediction of van der Waals complexes,[54−56] transition states,[57] and surface processes.[58,59] Accurate reference geometries and energies for the benchmark were generated by using the Cheap composite scheme (ChS)[60,61] and its recent jun-Cheap revision (jun-ChS),[16,32] with the latter appearing the best option because of the increased reliability for noncovalent interactions and the better description of the water dimer structure. Indeed, for (H2O)2, ChS and jun-ChS geometries were first compared to highly accurate CCSD(T)-F12b/CBS + fT + fQ + CV + REL + DBOC values.[62] The results, reported in Table S1 of the Supporting Information, show that bond lengths and valence angles are reproduced very accurately, with maximum errors of −0.003 Å and −0.2°, while there is a deviation of 3 °for the angle defining the orientation of the C2 axis of the acceptor water molecule with respect to the O–O axis. On the basis of the reliable geometry delivered by jun-ChS, this method was used as reference for both equilibrium geometries and electronic energies.

Preliminary B3LYP-D3/aug-cc-pVTZ computations of the HCN ⇌ HNC reactive PES were refined at the jun-ChS level. The nature of the identified stationary points (minima or saddle points) was checked through frequency calculations performed at each level of theory. All calculations have been carried out with the Gaussian software,[63] except the geometry optimizations at the ChS and jun-ChS levels, which have been performed using the CFOUR package.[64,65] Because revDSD-PBEP86 is not among the Gaussian built-in functionals, it has been defined by setting proper IOP flags on top of the DSD-PBEP86 functional.

Full geometry optimizations were performed for the complexes containing 2–4 H2O molecules, whereas for the 20 water model cut from the ice XI (010) surface, 8 molecules belonging to the cluster edge (see Figure S2 of the Supporting Information) were kept frozen at their positions in the crystal in order to prevent geometrical distortions causing a nonphysical breakdown of the crystalline pattern. The best-performing methods were employed within a QM/QM′ strategy for simulating the HCN ⇌ HNC isomerization on this cluster in order to evaluate the catalytic effect of the ice surface. For the purpose, we employed the ONIOM method[66] treating the reaction center (i.e., the adsorbate and four water molecules) at a higher level of theory (i.e., a double-hybrid DF or even jun-ChS), whereas a less computationally demanding method (i.e., a meta-hybrid DF) was used for the remaining molecules of the cluster. A much larger cluster containing 192 water molecules was also investigated by means of a three-layer (QM/QM′/MM) ONIOM approach enforcing the so-called mechanical embedding and employing the Amber force field.[67] In this case, the structural degrees of freedom of the adsorbate and the first 20 H2O molecules were optimized while freezing the coordinates of the remaining 172 water molecules to those of the regular (010) surface of ice XI. Test computations with the more refined electrostatic embedding showed negligible differences on the relative energies.

Rate constants were computed solving the multiwell one-dimensional master equation using the chemically significant eigenvalues method.[68] Rate coefficients were determined using conventional TST within the rigid-rotor harmonic-oscillator approximation,[69] also incorporating tunneling and nonclassical reflection effects by means of the Eckart model.[70] The rates evaluated at different temperatures were fitted by a simple Arrhenius equation or by the three-parameter modified Arrhenius equation proposed by Kooij[71,72]where A, n, and Ea are the fitting parameters, R is the universal gas constant, and the limiting Arrhenius behavior is recovered when n = 0. All the kinetic computations were performed with the MESS code.[68]

Results and Discussion

As widely discussed in the Introduction, the reliable modeling of interstellar ices is an extremely complex task, requiring the assessment of DFT methods for geometry and energy predictions that offer the proper balance between accuracy and computational burden. The lack of systematic studies addressing this issue for solid-state astrochemical processes calls for a dedicated benchmark. While small-size clusters cannot be fully representative of an extended substrate, the interaction of small molecules with water ice surfaces is generally guided by hydrogen bonds between the polar functional groups of the molecule and the exposed H and O atoms of the ice surface, which are already present in the smallest cluster models. Therefore, while the thermochemistry computed by using clusters composed of a small number of H2O molecules is not representative of real icy-grain chemistry, the outcomes of the benchmark are safely transferable to larger clusters. In the following subsections we report the results of our benchmark study, concerning first geometries and then reaction and activation energies. Finally, to scale-up to a more realistic water ice model, we report a full characterization of the PES of the HCN ⇌ HNC isomerization on clusters composed by either 20 or 192 H2O molecules.

Geometry Snow-Board

The HCN → HNC isomerization is an endothermic process involving a high activation energy and the jun-ChS results are close to the current best estimates[73] for both the reaction (61.6 vs 63.8 kJ/mol) and activation (198.5 vs 201.1 kJ/mol) energy. The addition of two water molecules leads to the formation of a hydrogen-bonded van der Waals adduct featuring the interactions between the H atom of HCN and the oxygen of one water molecule and between the N atom and one hydrogen of the second water molecule. Then, the reaction proceeds through a transition state for the (H2O)2-mediated proton transfer reaching, in this way, a post-reactive complex in which carbon is engaged in a weak H-bond with a hydrogen of the first H2O molecule, while the H atom of HNC interacts with the oxygen of the second water molecule. The structures of all the stationary points ruling the reactive PES are sketched in Figure together with selected geometrical parameters obtained at the jun-ChS level.

Figure 1

Stationary points on the reactive PES of the HCN ⇌ HNC isomerization catalyzed by two water molecules. Representative bond lengths (Å) obtained at the jun-ChS level are reported.

The accuracy of the considered DFT model chemistries has been evaluated with respect to jun-ChS values and the overall mean absolute errors (MAEs) and mean absolute relative errors (REs) have been evaluated over all the bond lengths, valence, and dihedral angles of the species involved in the PES. The full list of data can be found in Table S2 and Figure S1 of the Supporting Information. As a rule of thumb (with some exceptions for dihedral angles), triple-ζ basis sets show smaller errors than the corresponding double-ζ ones, with the improvement being less pronounced along the jun-, jul-, and aug-series. In general, the tested hybrid and meta-hybrid DFs on the one side, and the double-hybrids on the other, give similar trends for the MAEs, with the notable exception of BHLYP-D3 in conjunction with the jul-cc-pVDZ basis set, that strongly overshoots and the ωB97X-D functional that shows larger deviations from the jun-ChS reference values, especially for valence and dihedral angles. In the case of the BHLYP-D3/jul-cc-pVDZ model, MAEs as large as 0.09 Å, 6°and 10°were observed for bond lengths, valence, and dihedral angles, respectively. These results are related to the inability of reproducing a tight structure for the post-reactive complex. Specifically, one H-bond in CNH···(H2O)2 (see Figure ) is broken and the product collapses into an open structure. All in all, it can be observed that the most promising (meta-)hybrid DFs are PW6B95-D3, BMK-D3, M06-2X, and MN15 coupled to triple-ζ basis sets (or, at least, the jul-cc-pVDZ one). Concerning the double-hybrid functionals, the best structural predictions are delivered by DSD-PBEP86-D3 and revDSD-PBEP86-D3 that show comparable accuracy. In order to have a clearer picture of the performance of the different model chemistries in the prediction of the geometries involved in the HCN ⇌ HNC isomerization assisted by two water molecules, Figure reports the overall REs of each method, evaluated by averaging the REs of the geometrical parameters of all the species on the reactive PES. Inspection of this figure reveals that, among the (meta-)hybrid DFs, the best results for double-ζ basis sets are delivered by PW6B95-D3 and BMK-D3. In particular, PW6B95-D3/jul-cc-pVDZ, BMK-D3/aug-cc-pVDZ, and PW6B95-D3/aug-cc-pVDZ score REs in the 0.60–0.74% range. The PW6B95-D3 and BMK-D3 DFs are the best performers also in conjunction with triple-ζ basis sets showing REs around 0.55%. Concerning the double-hybrid functionals, it is apparent that their use in conjunction with a double-ζ basis set does not justify the computational overload in comparison with hybrid functionals; however, both DSD-PBEP86-D3 and its revision predict improved geometries when employed in conjunction with triple-ζ basis sets, reaching a RE of only 0.4% for the jul-cc-pVTZ basis set. In passing, it is interesting to point out that, these functionals have also demonstrated to be excellent performers in predicting structural and spectroscopic properties of gas-phase molecules.[74,75]

Figure 2

Total REs (%) of the geometries of the species on the PES of the HCN ⇌ HNC isomerization assisted by two water molecules for the investigated DFT methods with respect to jun-ChS reference values. For each functional, the different basis sets are reported in the following order: jun-DZ, jul-DZ, aug-DZ, jun-TZ, jul-TZ, and aug-TZ.

Skiing on Adsorption, Reaction, and Activation Energies

The functional/basis set combinations with the optimal accuracy/cost trade-off for geometry predictions have been identified in the previous section. Reactivity studies require the calculation of accurate formation and activation energies for the subsequent kinetic analysis. For this reason, some of the DFT methods delivering the best geometrical predictions have been selected and their accuracy for computing adsorption, activation, and reaction (electronic) energies explored using again jun-ChS results as references. In a first step, the impact of the geometry on the energetics has been assessed, by evaluating jun-ChS electronic energies for the different DFT structures. In a second step, the formation energies stemming from full DFT computations (for both geometries and energies) have been analyzed.

Electronic energies obtained at the jun-ChS level on top of selected DFT geometries are reported in Table , while the corresponding error analysis is presented in Figure .

Table 1

jun-ChS Formation Energies (kJ/mol) with Respect to Isolated HCN and (H2O)2 for Each Species along the HCN/HNC···(H2O)2 Isomerization PES Evaluated on Top of DFT Geometries

level of theory for geometry	HCN···(H2O)2	TS	CNH···(H2O)2	HNC + (H2O)2
PW6B95-D3/jul-DZ	–33.38	99.37	15.23	62.31
BHLYP-D3/aug-DZ	–33.40	99.21	14.87	61.88
PW6B95-D3/aug-DZ	–33.44	99.40	15.17	62.28
BMK-D3/aug-DZ	–33.12	99.13	15.04	62.26
M06-2X/aug-DZ	–33.32	99.67	15.38	62.28
MN15/aug-DZ	–33.39	99.36	15.61	62.63
PW6B95-D3/jul-TZ	–33.49	99.15	14.73	61.89
BMK-D3/jul-TZ	–33.47	99.18	14.59	61.66
M06-2X/jul-TZ	–33.38	99.49	14.91	61.85
MN15/jul-TZ	–33.51	99.33	14.98	62.03
DSD-PBEP86-D3/jul-TZ	–33.52	99.35	15.06	62.22
revDSD-PBEP86-D3/jul-TZ	–33.54	99.31	15.02	62.22
jun-ChS	–33.42	99.26	15.03	62.27

Figure 3

Error analysis for jun-ChS formation energy (kJ/mol) obtained on top of DFT geometries in comparison with full (both energies and geometries) jun-ChS results. Each color corresponds to a DFT model chemistry and collects absolute errors for the formation energy of each species along the PES with respect to isolated reactants: (1) pre-reactive complex; (2) transition state; (3) post-reactive complex; (4) products; and (5) MAE over all of the steps along the PES.

It is quite apparent that the energetic results obtained employing geometries optimized with all the tested methods are in remarkable agreement with the jun-ChS reference, with deviations smaller than 0.6 kJ/mol, even though some of them provide an unbalanced description of the different elementary processes. For example, the MN15/aug-cc-pVDZ and BMK-D3/jul-cc-pVTZ models yield excellent predictions of both the interaction energy of hydrogen cyanide with (H2O)2 and the transition state energy, with errors around 0.05 and 0.1 kJ/mol, respectively; however, the computed HNC formation energy (at the BMK-D3 level) and its interaction energy with the water dimer (at the MN15 level) show significantly larger errors. Among the (meta-)hybrid functionals, the best and most consistent energetic description is given by PW6B95-D3 in conjunction with jul- or aug-cc-pVDZ basis sets, which reaches an overall MAE (evaluated by considering the relative electronic energies of all the stationary points ruling the PES) close to 0.1 kJ/mol and a maximum deviation of 0.2 kJ/mol.

Moving to the double-hybrid DFs, the DSD-PBEP86-D3 and revDSD-PBEP86-D3 models in conjunction with the jul-cc-pVTZ basis set yield excellent performances, scoring a MAE of about 0.06 kJ/mol and reproducing the formation energies of all the elementary steps with a maximum deviation of 0.12 kJ/mol for the pre-reactive complex at the revDSD-PBEP86-D3/jul-cc-pVTZ level.

The relative electronic energies of all the stationary points fully evaluated at different DFT levels (i.e., energies and geometries) are collected in Table and the MAEs from the jun-ChS computations are shown in Figure . In general terms, the results mirror those obtained for jun-ChS energies evaluated on top of DFT geometries, with the only difference being the much larger deviations, which now span the 5–29 kJ/mol range. Furthermore, the relative stability of CNH···(H2O)2 is always strongly underestimated (becoming even negative with BHLYP-D3, BMK-D3, and MN15 functionals) except at the PW6B95-D3 and, especially, DSD-PBEP86-D3 and revDSD-PBEP86-D3 levels in conjunction with the jul-cc-pVTZ basis set. All the (meta-)hybrid DFs show MAEs larger than 10 kJ/mol, with the exception of PW6B95-D3, which is the only functional that reaches a MAE around 6 kJ/mol in conjunction with the jul- and aug-cc-pVDZ basis set and of 3.8 kJ/mol employing the jul-cc-pVTZ basis. The DSD-PBEP86-D3 and revDSD-PBEP86-D3 functionals confirm their good performances in conjunction with the jul-cc-pVTZ basis set, with MAE around 3 kJ/mol and maximum deviations of 10.4 kJ/mol. Hence, the model chemistries with the optimal accuracy for structural parameters are also the best choices for thermochemistry.

Table 2

DFT Formation Energies (kJ/mol) with Respect to Isolated HCN and Water Dimer (H2O)2 for Each Species along the HCN···(H2O)2 Isomerization PES

level of theorya	HCN···(H2O)2	TS	CNH···(H2O)2	HNC + (H2O)2
PW6B95-D3/jul-DZ	–34.12	86.96	8.34	56.74
BHLYP-D3/aug-DZ	–38.42	89.54	–1.35	51.12
PW6B95-D3/aug-DZ	–33.99	88.45	7.70	56.03
BMK-D3/aug-DZ	–36.30	81.76	–6.60	42.58
M06-2X/aug-DZ	–37.46	70.46	1.14	52.78
MN15/aug-DZ	–36.67	76.94	–2.82	47.31
PW6B95-D3/jul-TZ	–32.67	94.10	11.31	56.81
BMK-D3/jul-TZ	–35.97	88.16	0.11	48.23
M06-2X/jul-TZ	–37.16	76.00	4.64	54.06
MN15/jul-TZ	–35.26	85.53	1.92	48.88
DSD-PBEP86-D3/jul-TZ	–35.08	88.86	14.86	64.93
revDSD-PBEP86-D3/jul-TZ	–33.97	93.50	16.21	64.85
jun-ChS	–33.42	99.26	15.03	62.27

For both energy and geometry.

Figure 4

Error analysis for DFT formation energies (kJ/mol) in comparison with jun-ChS values. Each color corresponds to a DFT model chemistry (used for both geometry and energy) and collects absolute errors for the formation energy of each species along the PES with respect to isolated reactants: (1) pre-reactive complex; (2) transition state; (3) post-reactive complex; (4) products; and (5) MAE over all of the steps along the PES.

These results confirm the conclusions of recent benchmarks about the quality of PW6B95-D3/jul-cc-pVDZ and DSD-PBEP86-D3/jul-cc-pVTZ models for geometries, vibrational frequencies, and other spectroscopic parameters.[74,75] Noted is that the core–valence correlation has not been included for double hybrid functionals because it was not taken into account in their original parametrization and its contribution is anyway within the expected error bar at least for molecular systems containing only hydrogen and second-row atoms (see Table S3 of the Supporting Information for CV contributions in jun-ChS results). Furthermore, some test computations performed with quadruple-ζ basis sets showed that complete basis set extrapolation has a negligible effect on all the trends discussed above. For example, the relative electronic energies of the stationary points obtained by using the DSD-PBEP86-D3 functional in conjunction with the aug-cc-pVQZ basis set (ΔE = −35.04, 89.51, 14.83, and 64.78 kJ/mol) differ from the counterparts obtained employing the jul-cc-pVTZ basis set by 0.65 kJ/mol at most (for the TS). Finally, although triple-ζ basis sets possibly deliver more robust results for hybrid functionals, this computational level will be used in the following only to describe the environmental effects in the framework of QM/QM′ computations where the increased computational cost with respect to double-ζ results is not justified, in our opinion, by the marginally improved robustness.

Scaling-Up toward Extended Systems: Best Performers at Work

The benchmark performed for both geometries and energies permits the identification of the best candidates for setting up a QM/QM′ ONIOM strategy for the study of the HCN ⇌ HNC isomerization on large clusters capable of providing a more realistic modeling of the icy-grain and of the molecule–surface interactions.

At first, a cluster composed by 20 water molecules (shown in Figures and S2 of the Supporting Information) has been used, in which the pattern of exposed water molecules is suitable for a H-relay mechanism mediated by four water molecules. It should be noted that in ref (36) a proton-relay mechanism mediated by three water molecules, in turn solvated by seven additional waters, was used. In the present work, the four H2O molecules involved in the hydrogen transfer and the adsorbed species have been considered as the reaction center of the process under study; hence, they constitute the higher-level QM portion of the system. Following the outcomes of the benchmark study, the DSD-PBEP86-D3 functional in conjunction with the jul-cc-pVTZ basis set has been used for the purpose, while the remaining part of the cluster, treated at a lower QM′ level, has been described by the PW6B95-D3 DF in conjunction with the jul-cc-pVDZ basis set. The energetic profile of the HCN ⇌ HNC isomerization occurring on the (H2O)20 cluster is reported in Figure where it is also compared with that for the (H2O)2-mediated process. Going from the process assisted by two waters to that assisted by four water molecules in the (H2O)20 cluster lowers the energy of all the species present in the reactive PES. The most remarkable effect is the reduction of the energy barrier ruling the isomerization when considering the 20 water cluster in place of just two water molecules involved in the simplest possible relay mechanism.

Figure 5

Potential energy profile for HCN ⇌ HNC isomerization mediated by the (H2O)20 cluster and the (H2O)2 dimer. Red lines refer to the HCN isomerization catalyzed by (H2O)2 and both geometries and ΔE have been computed at the DSD-PBEP86-D3/jul-cc-pVDZ level. Black lines refer to the ONIOM results for the reaction catalyzed by (H2O)20. The ball and stick representation is used for atoms of the highest QM level (DSD-PBEP86-D3/jul-cc-pVDZ), while the tube representation is used for the atoms belonging to the QM′ (PW6B95-D3/jul-cc-pVDZ) portion. ΔE corrected for ZPVE are reported in parenthesis with ZPVEs calculated at the same level of theory as the corresponding energies and geometries.

The dependence of the energy profile on the number of water molecules involved in the relay mechanism was already pointed out.[35,36] However, only few water molecules were considered and no attempt to simulate the effect of ice bulk has been reported beyond the PCM level, whose reliability is, however, questionable for hydrogen-bonding solids. For comparison, Table lists the relative energies (corrected for the zero point vibrational energies, ZPVEs of the elementary steps obtained by Koch et al.,[36] and in the present work) (further details are given in Table S5 of the Supporting Information). As it can be seen, the relative energy for HCN interacting with the water cluster is only marginally affected by the cluster size, but there is a huge effect on the activation barrier. While an overall fair agreement between the present results and those obtained in ref (36), can be noted, there is a difference of about 18 kJ/mol for the energy of the transition state. This can be explained by considering that Koch et al.[36] investigated the role of the crystalline environment by optimizing for the different stationary points the positions of seven water molecules around the HNC···(H2O)3 complex without any constraint related to the arrangement of water molecules in icy structures. The importance of the morphological pattern in ice is highlighted by the present results: indeed, using a (H2O)20 cluster with the same molecular arrangement as in ice XI rules out the possibility of a process catalyzed by two or three water molecules. Rather, the molecular arrangement at the surface permits a process assisted by four water molecules (see Figure ).

Table 3

Relative Ground-State Energies (kJ/mol) with Respect to HNC···(H2O) Post Reactive Complex and Comparison with the Results of ref[36]a

	total H2O (n)	relay H2O (nR)	TSb	HCN···(H2O)nb
ref (36)	2	2	74.1	–42.3
	3	3	43.9 (−30.2)	–41.4 (−0.9)
	10	3	13.8 (−60.3)	–41.4 (−0.9)
B3LYPc	2	2	73.8	–39.8
	3	3	44.1 (−29.7)	–41.7 (−1.9)
PW6B95-D3d	2	2	70.9	–43.6
	3	3	52.4 (−18.5)	–42.0 (−1.6)
	4	4	49.3 (−21.6)	–42.7 (−0.9)
	192	4	36.6e (−33.4)	–40.0e (−3.6)
DSD-PBEP86-D3f	2	2	68.1	–51.5
	3	3	48.3 (−19.8)	–49.6 (−1.9).
	4	4	46.1 (−22.0)	–47.8 (−3.7)
	20	4	32.3g (−35.8)	–43.7g (−7.8)
			32.1h (−36.0)	–41.7h (−9.8)
	192	4	32.4i (−35.7)	–45.7i (−6.8)
jun-ChS	2	2	78.3	–50.0
	3	3	58.5 (−19.8)j	–48.0 (−2.0)j
	20	4	44.1 (−34.2)k	–40.8 (−9.2)k

Both the total number of water molecules (n) and the number of water molecules directly involved in the relay mechanism (nR) are indicated. All values include ZPVEs.

In parentheses is the difference with respect to (H2O)2 results.

6-31+G(d,p) basis set as in ref (36).

jul-cc-pVDZ basis set.

QM/MM energies and ZPVEs. 20 water molecules treated at the PW6B95-D3 level, the remaining molecules described by the Amber force field.

jul-cc-pVTZ basis set.

ONIOM geometries and ZPVE. DSD-PBEP86/jul-cc-pVTZ for adsorbate and molecules involved in the relay mechanism, PW6B96-D3/jul-cc-pVDZ for the water molecules not involved in the relay mechanism.

DSD-PBEP86/jul-cc-pVTZ energies with geometries and ZPVE at the PW6B95-D3/jul-cc-pVDZ level.

DSD-PBEP86:PW6B95-D3:Amber energies on PW6B95-D3:Amber geometries. ZPEs at PW6B95-D3:Amber level.

jun-ChS electronic energy, PW6B95-D3/jul-cc-pVDZ geometry, and ZPVE.

jun-ChS:PW6B95 electronic energy, PW6B95-D3/jul-cc-pVDZ geometry, and ZPVE.

Further support to the reliability of the results is provided by the comparable barrier obtained by another ONIOM computation in which the high-level part of the system [HCN···(H2O)4] is treated at the jun-ChS instead of DSD-PBEP86-D3 level without any additional geometry optimization (last line of Table ). What is even more gratifying is that the differences between the results obtained for the smallest HCN···(H2O)2 model and the larger model clusters (values in parenthesis in Table ) obtained at the DSD-PBEP86 level are in quantitative agreement with the jun-ChS counterparts. This paves the route toward the computation of very reliable parameters for reactions occurring on icy grains by combining jun-ChS results for small models and ONIOM(DSD-PBEP86:PW6B95-D3) values for large model clusters.

This approach can be further extended to very large models by employing a three-layer QM/QM′/MM ONIOM model. In order to also analyze this aspect, we have embedded the HCN···(H2O)20 cluster in a large model of ice-XI containing 172 water molecules described by the Amber force field (see Figure ). The results collected in Table show that inclusion of the MM layer further stabilizes the HCN isomer with respect to the HNC counterpart by about 4 kJ/mol, but has a negligible effect on the energy barrier (less than 0.4 kJ/mol). Taking into account the estimated error bar of the overall computational approach (about 4 kJ/mol), the results obtained for the HCN···(H2O)20 model can be considered essentially converged with respect to further extension of the ice substrate.

Figure 6

Structural model for the (H2O)192 cluster treated by three-layer ONIOM DSD-PBEP86:PW6B95:Amber strategy (geometry at PW6B95:Amber level). Ball and stick and tubular representation for the QM sections treated at the DSD-PBEP86/jul-cc-pVTZ and PW6B95-D3/jul-cc-pVDZ levels, respectively.

Reaction Rates

In an astrochemical context, HNC can either isomerize to HCN or diffuse on ice surfaces and then react with another molecule (e.g., CH2NH to produce acetonitrile) at the low temperatures typical of the ISM.

The reaction rates computed for the HNC ⇌ HCN isomerization with the methodology described in Section are shown in Figure . It is apparent that the rates computed for the HNC···(H2O)2 system [Figure , panels (c,d)] are very small irrespective of the inclusion or not of tunneling. The situation is completely different for the HNC···(H2O)20 model, where the rate not including tunneling (corresponding to the one used by Koch and coworkers[36]) remains very small at low temperatures (Figure b), but inclusion of tunneling (Figure a) permits an effective reaction even at temperatures characteristic of the ISM. Noted is that the rates computed taking tunneling into account show a clear bimodal shape and cannot be fitted by a simple Arrhenius (or Kooij) function.[71,72]

Figure 7

Reaction rates for the HNC ⇌ HCN isomerization including (Eckart) or excluding (no tun) tunneling. Panels (a,b) refer to the HNC···(H2O)20 model, whereas panels (c,d) refer to the HNC···(H2O)2 model.

Unfortunately, the diffusion coefficients of HNC (or even HCN) on ice have not yet been reported.[76] According to a recent classification of ice adsorbates,[77] HCN (hence probably HNC) is assigned to the intermediate class, which induces some deformation of the surface, but does not form hydrates nor penetrates rapidly into the ice bulk. An upper limit to the surface diffusion coefficient can be estimated with reference to the guess of 4 × 10–11 cm2 s–1 at 130 K provided by Livingston et al. for SO2,[78] which corresponds to a mean distance of 1260 Å in 1 s. Because the computed isomerization rate at 130 K is about 10–4 s–1 (which lowers to 1 × 10–10 s–1 at 50 K), the average diffusion of HNC before isomerization can reach 100 Å at 130 K (1 cm at 50 K).

Therefore, if the formation of aminoacetonitrile is faster than the isomerization to HCN when HNC and CH2NH are nearest neighbors,[34] our results suggest that diffusion of HNC along significant distances could permit the formation of aminoacetonitrile on icy grains containing CH2NH even at low concentrations.

Conclusions and Outlook

The main aim of this work was the implementation and validation of a general computational strategy for the study of the thermochemistry and kinetics of chemical processes taking place on interstellar icy-grains. To this end, composite methods rooted in the coupled cluster ansatz have been combined with hybrid and double hybrid functionals together with molecular mechanics force field to characterize the stationary points ruling the reactive potential energy surfaces on model clusters sufficiently large to minimize spurious boundary effects. Next powerful master equation/TST models have been employed to compute reaction rates including tunneling effects. As a demanding test case, we have selected the HCN/HNC reactions for which the available computational results are not fully satisfactory.

Ten different (meta-)hybrid and double-hybrid density functionals have been considered in conjunction with the jun-, jul-, and aug-cc-pVnZ basis sets of double- and triple-ζ quality, and their accuracy in predicting geometries together with thermochemical and kinetic data (adsorption, activation, and reaction energies) has been assessed in comparison to reference values computed using the jun-ChS composite method. This benchmark has led to the conclusion that, among (meta-)hybrid functionals, BMK-D3 and PW6B95-D3 in conjunction with partially augmented double- and triple-ζ basis sets yield the most reliable description of geometries, with the optimal trade-off between accuracy and computational cost being offered by the PW6B95-D3/jul-cc-pVDZ model chemistry. Concerning double-hybrids, DSD-PBEP86-D3 and revDSD-PBEP86-D3 in conjunction with the jul-cc-pVTZ basis set deliver accurate predictions of both geometries and reaction energies. Next, these outcomes have been used to investigate the effect of cluster size and ice surface on the isomerization process of HCN. In particular, a cluster containing 20 water molecules has been cut from the (010) surface of ice XI and used in a multiscale ONIOM calculation, in which the reaction center has been modeled at the DSD-PBEP86-D3/jul-cc-pVTZ level, while for the remaining portion of the (H2O)20 cluster, the PW6B95-D3 functional has been employed in conjunction with the jul-cc-pVDZ basis set. This approach has allowed the proper modeling of the surface with an accurate yet cost-effective strategy. The pivotal role of the structural arrangement of surface molecules in driving the evolution of catalytic processes has been pointed out. The accuracy of the results has been further improved by combining jun-ChS results for small models to QM/QM′ (DSD-PBEP86:PW6B95-D3) values for medium-size model clusters and/or three-layer QM/QM′/MM computations for very large clusters.

On top of these computations, reaction rates have been computed by methods rooted in the transition state theory including tunneling which plays the dominant role at low temperature for processes involving the motion of light atoms. At variance with previous investigations, our results show that the isomerization is ruled by a proton relay mechanism directly involving four water molecules, but tuned by relatively distant waters belonging to the model cluster employed to mimic the ice surface. The resulting activation energy is strongly reduced with respect to that governing the isomerization of the bare HCN molecule, but only tunneling allows for effective isomerization of HNC in the harsh conditions characterizing astrochemical processes.

Together with the intrinsic interest of the studied system, the results of the present work have allowed to define the best strategy for future modeling of iCOMs-ices interactions in the framework of a QM/QM′/MM approach. This also represents the starting point for hybrid QM/QM′/periodic approaches, in which the outcome of the multiscale (QM/QM′) description of the cluster is corrected for environmental effects obtained by simulating the surface using periodic boundary conditions.[79] However, the crystalline water ice surfaces usually employed to simulate icy dust grains could be inadequate to describe their amorphous structure. Work in this and related directions is under way in our laboratory in order to achieve a more realistic modeling of chemical processes occurring on icy mantles of interstellar grains.