QM/MM Study of the H2 Formation on the Surface of a Water Ice Grain Doped With Formaldehyde: Molecular Dynamics and Reaction Kinetics

Formaldehyde has been widely observed in the icy mantle of interstellar grains. H2CO may be formed from successive hydrogenations of CO and may further contribute to the chemical complexity of the Interstellar medium (ISM) participating to heterogeneous reactions with colliding gas phase atoms. Within this context, Eley-Rideal and Langmuir-Hinshelwood rate constants of H2 formation on a formaldehyde doped amorphous water ice grain model of the ISM, were computed over a wide temperature range [15–2000 K]. We used classical molecular dynamics (MD) simulations to build the model of the H2CO doped ice surface. Then we studied theoretically by means of hybrid QM/MM ab initio and molecular mechanics methodology (ONIOM) H atoms abstraction from formaldehyde molecules and the H2 formation. Specifically, we investigate the reactivity of the gas phase H atom toward one formaldehyde molecule lying at one of the slab surfaces. The reaction path and the energetics are predicted, the mechanism is found to be exothermic by 14.89 kcal/mol and the barrier is 6.75 kcal/mol at the QM level CBS/DLPNO-CCSD(T)//ONIOM/aug-cc-pVTZ. We employ two approaches that take into account tunnelling and non-classical reflection effects by means of the Zero Curvature Tunnelling (ZCT), and the Small Curvature Tunnelling (SCT) which all provided comparable results to predict the kinetics of the reaction path. The rate constants show important quantum tunnelling effects at low temperatures when compared to rates obtained from the purely classical transition-state theory (TST) and from the canonical variational transition state theory (CVT). Corner cutting effects are highlighted in the SCT calculations by 4 to 5 orders of magnitude with respect to ZCT rate constants at low temperatures.

Formaldehyde has been widely observed in the icy mantle of interstellar grains. H 2 CO may be formed from successive hydrogenations of CO and may further contribute to the chemical complexity of the Interstellar medium (ISM) participating to heterogeneous reactions with colliding gas phase atoms. Within this context, Eley-Rideal and Langmuir-Hinshelwood rate constants of H 2 formation on a formaldehyde doped amorphous water ice grain model of the ISM, were computed over a wide temperature range [15-2000 K]. We used classical molecular dynamics (MD) simulations to build the model of the H 2 CO doped ice surface. Then we studied theoretically by means of hybrid QM/MM ab initio and molecular mechanics methodology (ONIOM) H atoms abstraction from formaldehyde molecules and the H 2 formation. Specifically, we investigate the reactivity of the gas phase H atom toward one formaldehyde molecule lying at one of the slab surfaces. The reaction path and the energetics are predicted, the mechanism is found to be exothermic by 14.89 kcal/mol and the barrier is 6.75 kcal/mol at the QM level CBS/DLPNO-CCSD(T)//ONIOM/aug-cc-pVTZ. We employ two approaches that take into account tunnelling and non-classical reflection effects by means of the Zero Curvature Tunnelling (ZCT), and the Small Curvature Tunnelling (SCT) which all provided comparable results to predict the kinetics of the reaction path. The rate constants show important quantum tunnelling effects at low temperatures when compared to rates obtained from the purely classical transition-state theory (TST) and from the canonical variational transition state theory (CVT). Corner cutting effects are highlighted in the SCT calculations by 4 to 5 orders of magnitude with respect to ZCT rate constants at low temperatures.

INTRODUCTION
H 2 molecules are ubiquitous in the Universe and constitute the most abundant species and consequently are the main trigger of dust chemistry which is responsible for the formation of key molecules, such as H 2 O, NH 3 , and CH 3 OH. Detection of H 2 was achieved in the interstellar medium (ISM) via a rocket flight in 1970 (Carruthers, 1970). It has been shown that almost certainly H 2 forms at the surface of dust grains and icy mantles in a variety of interstellar environments (Hollenbach and Salpeter, 1970;Williams et al., 2003;Vidali, 2013;Wakelam et al., 2017). Extensive computational (Katz et al., 1999;Cuppen et al., 2010;Cazaux et al., 2011), and experimental studies (Pirronello et al., 1997;Creighan et al., 2006;Watanabe et al., 2010) over the 20 past years have been conducted. Many focused on testing the role of grain models on the energetics prediction of the chemical routes leading to H 2 formation.
Among the grain models, there appear carbonaceous (Farebrother et al., 2000;Ivanovskaya et al., 2010), and silicate ones (Goumans, 2011;Kerkeni and Bromley, 2013). From the modeling point of view molecular clusters and also periodic slabs have served to mimic interstellar grains (Farebrother et al., 2000;Ivanovskaya et al., 2010;Goumans, 2011;Chaabouni et al., 2012;Navarro-Ruiz et al., 2014;Oueslati et al., 2014;Navarro-Ruiz et al., 2015;Oueslati et al., 2015;Navarro-Ruiz et al., 2016;Song and Kästner, 2017), upon which H 2 formation was investigated. Several experimental techniques e.g., (Latimer et al., 2008;Wakelam et al., 2017) at cryogenic temperatures have been realized to investigate the efficiency of H 2 formation in the laboratory. Given the very low temperatures ≈10 K in dense regions of space, dust grains are covered with layers of icy mantles (Ehrenfreund and Schutte, 2000). These serve as a reservoir for molecular accretion, that would subsequently allow the reagents to thermalize, diffuse and react (Congiu et al., 2014), and finally desorb in a later phase during star formation processes. These formed molecules are detected in warmer diffuse clouds of the interstellar medium. The essential enablers of this rich and complex reactivity at the surface of grains may rely in some cases on radical-radical recombination or the presence of UV irradiation and cosmic rays. Dense molecular clouds regions are however less exposed to radiations, suggesting that the grain may truly act as a catalyst reducing the reaction barriers.
Pure H 2 O ice was detected at 3 μm by Gillett andForrest since 1973 (Gillett et al., 1973), it is the most abundant molecular ice in dense ISM regions. While formaldehyde absorption features in ices were observed in young stellar objects (YSOs) in the 5.7-6 μm region. These detections hint to abundances around 6% with respect to H 2 O (Gibb et al., 2004;Boogert et al., 2008).
Formaldehyde (CH 2 O) ices are believed to form upon water-CO mixtures hydrogenation (Hiraoka et al., 1994;Hiraoka et al., 2002;Watanabe and Kouchi, 2002;Fuchs et al., 2009;Pirim et al., 2010). Comets are also known to contain water related ices in their cores, and when heated by nearby sources, can participate to a complex chemistry network (Schutte et al., 1993;Bossa et al., 2009;Vinogradoff et al., 2011). Laboratory experiments have tackled thermal evolution of comets in order to understand chemical pathways towards Complex Organic Molecules (COMs) formation (Duvernay et al., 2014).
Despite the past extensive investigations, some uncertainties remain, such as the exact composition of the grain models, and the role of mixed ices in the efficiency of H 2 formation. The chemical mechanisms of H 2 formation and the underlying energetic pathways can be model-dependent, and can in some cases present hindering energetic barriers preventing H 2 formation at the low temperatures of the ISM. It is thus timely to provide a detailed energetic description of the possible reaction mechanisms on each specific grain model. In order to be used by astrochemical models, diffusion barriers, adsorption energies and analytic expressions of rate constants should be predicted.
Hydrogen addition and abstraction from formaldehyde have been investigated computationally with DFT by Goumans (Goumans, 2011), and CCSD (T) (by one of the authors) (Siaï et al., 2016) rate constants predictions in gas phase, and on pure water ice (Song and Kästner, 2017). In Woon (2002), quantum chemical electronic structure calculations of H + H 2 CO reaction in gas phase and in the presence of water clusters were lately investigated, the catalytic activity of the ice models was analyzed as a function of the activation energies magnitude and without any kinetic modeling. Song and Kästner (2017), employing a QM/ MM model, studied the H 2 formation on an amorphous solid water (ASW) upon which one formaldehyde molecule was adsorbed, and rate constants for both Langmuir-Hinshelwood and Eley-Rideal mechanisms were evaluated with the instanton theory, over the temperature range 60-300 K. On the other hand, Monte Carlo simulations (Cuppen et al., 2009) of formaldehyde formation upon CO hydrogenation on interstellar ice surfaces show that H 2 CO molecules are found to predominantly exist in the deeper layers of their ice model. H abstraction and addition reactions were also studied experimentally (Chuang et al., 2016) in mixed CO, CH 3 OH, and H 2 CO ices, where the authors confirm the occurrence of the product channels leading to CH 3 OH through hydrogenation, however successive abstraction reactions lead to the formation of HCO and back to CO molecules.
This study may serve as a useful investigation for further model variants to be explored computationally, and for suggesting possible constraints on new experimental setups relevant to ice models and H 2 formation on icy grains in general.
In the present paper we are specifically interested in studying H 2 formation from a formaldehyde-doped ice sample prepared by Molecular Dynamics techniques. The H abstraction from ice formaldehyde by a free incoming gaseous H atom will be investigated considering both reaction mechanisms: -The Eley-Rideal (ER) mechanism: -The Langmuir-Hinshelwood (LH) mechanism: Where H (gas), H 2 (gas) are the free hydrogen atom and molecule respectively in gas phase, while H (ads), H 2 (ads) are the adsorbed species on the grain model. We use a slab model of formaldehyde doped water ice, and choose a reactive site where an accessible CH 2 O molecule can be reached upon reaction. We investigate the reaction energetics for H 2 (gas) formation with QM/MM methodologies. Semi-classical rate constants at a large temperature range are predicted using sophisticated Variational Transition State Theory (VTST) (Fernandez-Ramos et al., 2007), that allow unlike other methods to compute rate constants for temperatures above the cross over temperature T c and also at very low temperatures ≈15 K characteristic of dense regions in the ISM where molecules are believed to form. We compare our results to available data from the literature and investigate the role of the present model grain and computational methodologies on the prediction of the energetics and the rate constants. In Computational Details we describe the methodology to generate the doped ice slab, the adopted computational methods for electronic structure and rate constants calculations. Results and their implications are discussed in Results and Discussion. Finally, the main results are put into perspective in Conclusion.

Classical Molecular Dynamics Calculations
The mixed formaldehyde/water ice was generated by classical molecular dynamics calculations using the Gromacs code (Abraham et al., 2015). The water molecules were described with the TIP4P/2005 model that proved to reproduce fairly well ice along the phase diagram (Abascal and Vega, 2005;Michoulier et al., 2018). OPLS force field (Jorgensen et al., 1983;Cornell et al., 1995) was employed for the formaldehyde molecules. The combination of OPLS and water models derived from TIP4P has been validated in Refs. (Collignon and Picaud, 2004;Ballenegger et al., 2006). The equations of motion were integrated by the leapfrog algorithm with a 1 fs time step. The Lennard-Jones and the short-range part of the Coulomb interactions were truncated at 10 Å. The long-range part of the Coulomb interaction was evaluated using the smooth particle-mesh Ewald method with a relative tolerance of 10 -5 , fourth-order cubic interpolation and a Fourier spacing parameter of 0.12. Water and formaldehyde molecules were treated as rigid bodies using the SETTLE algorithm.
Amorphous ice was doped with formaldehyde at the ratio 20 H 2 CO:1000 H 2 O (Pontoppidan et al., 2004). The slab is formed from 20 CH 2 O randomly distributed in a box containing 1000 H 2 O molecules. In order to favor diffusion of the molecules, an annealing procedure is carried out in the NPT ensemble using a Berendsen thermostat, heating the sample from 100 to 300 K and then cooling it back to 110 K for a total duration of 10 ns. The obtained box, of dimensions 26 × 30 × 35 Å, is then elongated along the z direction to form two interfaces exposed to vacuum. The resulting box cell dimensions are 26 × 30 × 100 Å. A final equilibration in the NVT ensemble was then performed at 110 K during 10 ns. A side view of the doped ice sample is shown in Figure 1A. The amorphous character of water was insured by checking the water density as shown in Figure 1B. Moreover, the compacity and low porosity of the ice are not affected by the presence of the organic molecules, the density of the ice slab is only changing from 0.955 g/cm 3 for pure amorphous ice to 0.989 g/cm 3 for our doped ice sample. H 2 CO molecules are distributed all over the slab but few of them remain at the surface. The abstraction reaction is subsequently studied from one of these surface formaldehyde molecules within a QM/MM scheme. It should be noted that at the low concentration considered here, the mutual interactions between H 2 CO molecules is screened by water molecules. The water network is also only slightly affected by the presence of the H 2 CO molecules.

Quantum Chemical Calculations
Because of its amorphous nature and large molecular size, the whole system cannot be studied within a quantum mechanical treatment. Therefore, it has been divided into two parts. The first one, where the electronic energy of the system is calculated with an ab initio framework (QM), is located around the site where the reaction occurs. On the other hand, classical mechanics methods (MM) were used for the rest of the icy layer. We employ the hybrid QM/MM ONIOM methodology (Chung et al., 2015) where the active species involved in the reaction process, i.e., H, CH 2 O and 10 H 2 O within a 4 Å radius of the CH 2 O molecule are studied with the hybrid meta-GGA functional M05-2X (Zhao et al., 2006) combined with an augmented correlation-consistent polarized valence basis set of Dunning namely, the aug-cc-pVTZ (Dunning, 1989;Kendall et al., 1992). This approach satisfactorily predicts the geometry and is particularly suitable for treating weak interactions as well as for providing accurate activation energy barriers (Hohenstein et al., 2008;Leverentz et al., 2013;Sirianni et al., 2018). M05-2X has been tested to perform in a similar way as M06-2X in the case of hydroxylated pyroxene clusters (Kerkeni et al., 2019). The rest of the system, i.e., the icy grain minus CH 2 O + 10H 2 O is treated with the Amber MM force field. The mechanical embedding scheme (Chung et al., 2015) was employed for all calculations. For open-shell structures, unrestricted wave functions were employed. All the DFT quantum chemical calculations were performed with the Gaussian 16 code (Frisch et al., 2016).
Initially the CH 2 O + 10H 2 O molecules in the QM region were allowed to relax in a QM/MM calculation, the final structure was afterwards used while freezing the 10 H 2 O molecules. Therefore, the reactive H and CH 2 O fragments were allowed to relax without any constraint during the reaction mechanism. The remainder of the grain model was kept frozen during all calculations. Geometry optimizations were carried out using the default Gaussian convergence thresholds. The nature of the stationary points of the potential energy surface (PES), minimum or saddle point, was checked by computing unscaled harmonic frequencies.
Transition states (TSs) were searched for using the Berny algorithm (Schlegel, 1982). The Minimum Energy Path (MEP) was computed using the Intrinsic Reaction Coordinate (IRC) method (Fukui, 1981;Hratchian and Schlegel, 2005) with a gradient step size of 0.01 a 0 .(amu) 1/2 . For 46 points along each side of the MEP, the harmonic frequencies were computed at the M05-2X/aug-cc-pVTZ/Amber level. The MEP energies were finally corrected by adding the unscaled zero-point vibrational energy (ZPE), leading to the vibrational adiabatic potential energy V G a . In order to obtain more accurate energies at an affordable cost, the DFT results were improved using a methodology recently used to obtain accurate binding energies of atoms and small molecules at the surface of both crystalline and amorphous ices (Duflot et al., 2021). To this end, the optimized geometries of the QM zone for all stationary and IRC points were recovered. The resulting cartesian coordinates were then used as inputs for a CBS/DLPNO-CCSD(T) (Guo et al., 2018) single point calculation, using the ORCA code (Neese, 2018;Neese et al., 2020). Resolution of the identity (RIJK) was used (Weigend et al., 2009) with the auxiliary basis set from Weigend (Weigend, 2008). For DLPNO-CCSD(T), the normal criteria were used as discussed by Liakos et al. (2015). The frozen core approximation was also used so carbon and oxygen 1 s orbitals were not included. The auxiliary basis set for correlation is taken from Hellweg et al. (2007). For Complete Basis Set (CBS) extrapolation, def2-TZVPP and def2-QVZPP basis sets were employed (Weigend and Ahlrichs, 2005). The Hartree-Fock and the correlation energies were extrapolated separately as shown in (Neese and Valeev, 2011) with α = 7.880 for Hartree-Fock and β = 2.970 for correlation. Interestingly, the T 1 diagnostic (Lee et al., 1995) remains small along the path, with a maximum value of 1.2% at the TS geometry, thus confirming the absence of multireference character to the wavefunction. This diagnostic was slightly larger for gas phase species, peaking at 2.4% for the TS. By replacing the energy part of the total ONIOM calculation by the CBS/ DLPNO-CCSD (T) energies, a new ONIOM extrapolated energy is obtained. On the other hand, the DFT frequencies and the resulting ZPE values were kept unmodified. This procedure was applied on the stationary and all the IRC points along the calculated reaction path. Frontiers in Astronomy and Space Sciences | www.frontiersin.org March 2022 | Volume 9 | Article 807649

Conventional Transition State Theory
The conventional transition-state-theory rate constant is the product of a symmetry number by k B T/h multiplied by the ratio of a partition function evaluated at the saddle point to one evaluated at the reactant, times an exponential decay of the classical barrier height as in Eq. 3. The partition functions are evaluated as products of electronic and vibrational partition functions. As we are dealing with surface reactions, we discard the relative translational partition function per unit volume and the rotational ones. The conventional transition state-theory rate constant, k TST (T), is given by the following expression (Truhlar et al., 1985): where ΔV a is the forward classical barrier height, h is the Planck constant, and k B is the Boltzmann constant. Q ‡ and Q R are the TS and the reactant partition functions, respectively, which are relevant to internal vibrational motions, and the electronic distribution. All vibrational modes are treated as separable harmonic oscillators.

Canonical Variational Transition State Theory
Canonical variational transition-state theory (CVT) (Liu et al., 1993) has been used to estimate the rate constants. The canonical variational, transition state theory rate constant k CVT is obtained by finding the location of the minimum flux through the transition state at a given temperature. Note that in using a free energy of activation criterion to determine the location of the variational transition state, both entropic and energetic contributions are taken into account. We may write: where V MEP (s) is the classical potential along the minimum path (MEP) referenced to the reactant energy, s is the mass scaled reaction coordinate, and s* ,CVT is the location of the canonical transition state on the reaction path, i.e., where the flux is minimal. s = 0 corresponds to the transition state; at the reactant s is a finite negative number, while for the product s is a finite positive number.

Semiclassical Reaction Probabilities
Light atoms such as hydrogen tunnel through barriers that cannot be surmounted classically. This quantum effect can be important at low temperatures. The exponential damping of wave tunnelling through a given barrier leads to the tunnelling probability P(E) at total energy E: where θ(E) represents the amount of exponential decay. This function may be obtained in several ways.

Zero Curvature Tunnelling Method
When the reaction path curvature can be neglected, the proper tunnelling path is the MEP. For tunnelling along the MEP, θ(E) is the imaginary action integral given by: V AG is the maximum of the ground state vibrational adiabatic potential energy: where V G a (s) V MEP (s) + ZPE(s), V MEP is the classical potential energy of the MEP and ZPE is the zero-point energy correction relevant to the frequencies orthogonal to the reaction path. The s 0 and s 1 are the reaction coordinates of the classical turning points in the reactant and product valleys, respectively, TABLE 1 | Calculated exothermicity and barrier height compared to gas phase results (kcal/mol). The inclusion of ZPE is denoted by a star symbol.

Grain
Gas and μ is the reduced mass. This approach is called the zero curvature tunnelling (ZCT) method (Skodje et al., 1981).
In the ZCT approach, tunnelling is calculated from the semiclassical (i.e., WKB) ground-state transmission coefficient κ(T) (Skodje et al., 1981): The transmission coefficient κ is temperature dependent, which takes into account tunnelling and non-classical reflection effects, κ(T) is given by Small Curvature Tunnelling Method Liu et al. (1993) have reported that tunnelling probabilities may increase because of the corner-cutting effect, which is not accounted for in the ZCT method. More advanced tunnelling methods are based on the inclusion of deviations between the tunnelling path and the reaction path (Vandeputte et al., 2010). In the present work, the centrifugal dominant small curvature semiclassical adiabatic ground state (CD-SCSAG) tunnelling method (Lu et al., 1992;Liu et al., 1993) , is applied, where the effect of the reaction curvature is included by replacing the reduced mass μ by an effective mass μ eff (Liu et al., 1993). The tunnelling probability at energy E is: where θ(E) is the imaginary action integral along the tunnelling path as in Eq. 6. κ(T) is calculated as in Eq. 9 and the thermal rate constant is tunnelling corrected in the corresponding way as done for ZCT. The TST, CVT, CVT/ZCT and the CVT/SCT rate constants are computed using the POLYRATE program, version 2010-A . The CVT rates were fitted using the analytical expression proposed by Zheng and Truhlar (2010): where α, β, γ and T 0 are constants optimized to match the calculated rates.

Gas Phase Abstraction
In order to validate the methods employed in the ONIOM calculations, namely M05-2X/aug-cc-pVTZ and CBS/ DLPNO-CCSD(T), we have employed them to calculate the gas phase reactions. It should be kept in mind that the CBS single point calculations were performed at the M05-2X geometries and ZPE is evaluated at this DFT level too. The results are compared to previous works (Goumans, 2011;Ruscic, 2015;Siaï et al., 2016;Song and Kästner, 2017) in Table 1. The computed energetics given in this table are the exothermicity (ΔE) and the vibrational adiabatic barrier (ΔV a ) between the transition state and the reactants and between the products and the reactants, respectively. The corresponding ZPE-corrected values are labelled by a star symbol. It can be seen that the M05-2X/aug-cc-pVTZ overestimates the barrier height and underestimates the exothermicity when compared to the CBS/DLPNO-CCSD(T) values. The same behavior was also found previously at the MP2/cc-pVTZ level (Siaï et al., 2016). The CBS/DLPNO-CCSD(T) gas-phase barrier height (5.85 kcal/mol) is slightly below the previous CCSD(T)/cc-pVTZ//MP2/cc-pVTZ result of 6.69 kcal/mol (Siaï et al., 2016) and the CCSD(T)-F12a/cc-pVTZ-F12 value of 6.68 kcal/mol. The non-ZPE-corrected barriers can also be compared to the PWB6K-D3/def2-TZVP and (U)CCSD(T)-F12/def2-TZVP results of Song and Kästner (2017), i.e., 6.60 and 7.03 kcal/mol, respectively. As shown by these authors  [see Table 4 in (Song and Kästner, 2017)], the barrier height is very sensitive to the choice of the functional and basis set, leading sometimes to negative values, but always below the (U)CCSD(T)-F12 result. Hence, our M05-2X/aug-cc-pVTZ value of 9.08 is too large. On the other hand, the present CBS/ DLPNO-CCSD(T) value of 7.62 kcal/mol is only about 8% larger than the (U)CCSD(T)-F12 one. These results confirm that the present improved treatment, i.e., the composite methodology beyond DFT, ensures reliable results on the energetics of the reaction.

Stationary Points
We initiated the reaction process by the reactive collision of one single gaseous H atom with an already adsorbed CH 2 O molecule dangling at the doped ice grain surface (Figure 2). The present geometry is similar to the most abundant one obtained by Song and Kästner on pure ASW ice (see Figure 1A in ref. (Song and Kästner, 2017). It is worth mentioning that in ref (Song and Kästner, 2017) a single formaldehyde molecule was adsorbed on their ice models. In the present model, the reactive formaldehyde oxygen atom forms a hydrogen bond with the nearest water hydrogen atom with a length of about 2.00 Å, close to the 1.90 Å value obtained by Song and Kästner. A value of 1.88 Å was also obtained at BHLYP/6-311++G (d,p) level on a 18-water cluster by Rimola et al. (Rimola et al., 2012). The search of the transition state corresponds to the geometrical configuration where the hydrogen in CH 2 O is abstracted by the incoming H atom ( Figure 3 and Table 2). At the transition state the CO bond is shortened by ≈0.017 Å with respect to the reactant one (see Figure 3 and Table 2).  (Siaï et al., 2016) and on pure ASW ice (Song and Kästner, 2017). The relevant geometrical parameters are given in Figure 4 and Table 2 for all involved species in the reaction mechanism. The calculated H-H distance at the TS (1.068 Å) is slightly larger than the 1.022 Å obtained by Song and Kästner (2017). All attempts to locate reactive and product Van der Waals complexes failed. Song and Kästner (2017) mention a complex with an energy very close to the reactants but no details are given about its geometry. In the gas phase case (Siaï et al., 2016), it was concluded that the product complex was electronically unstable after including BSSE effects. The reactant complex, with the incoming H atom at about 3.35 Å, lies very close to the reactant (-023 kcal/mol at CCSD(T) level) and did not play any significant role in the dynamics study of the abstraction reaction. Figure 5 shows a comparison between the geometries of the stationary points and the last points of the IRC calculations. Beyond these geometries, the IRC calculation stops due to the stability of the energy. On the reactant side, at CBS level, the last IRC point lies 0.23 kcal/mol (with ZPE) above the energies of the isolated grain + H system. This suggests that the reaction barrier for the Eley-Rideal mechanism would be 0.23 kcal/mol higher than that for the Langmuir-Hinshelwood one (see Table 4). The H atoms is at ≈ 3.4 Å from the CH 2 closest hydrogen. On the other side, the last IRC point is 1.71 kcal/mol above the isolated products (see also Figure 6). The closest H atom of H 2 is ≈2.8 Å from the carbon atom. For both cases there remains a large difference between the orientation of the isolated H 2 CO ( Figure 5A) or HCO ( Figure 5B) molecules compared to the last IRC point, which proves that the convergence of the MEP 3 | Calculated harmonic vibrational frequencies (in cm −1 ) and zero-point energies (in Hartrees) for the stationary points.

Species
Frequencies ZPE H 2 CO 42, 53, 68, 91, 151, 198, 1246, 1285, 1553, 1845, 3015, 3106 0.028822 TS i1737, 40, 46, 48, 94, 99, 136, 341, 440, 1184, 1223, 1322, 1405, 1891, 3024 0.025724 HCO 36, 52, 62, 82, 137, 298, 1132, 1968, 2876 0.015134 H 2 4533 0.010327 calculation has already reached local minima in both sides. The same behavior is also found for the TS as shown in Figures 5C &  5D, where one can see the orientation of the reactant and product compared to the TS. Several attempts were made to determine whether isolated reactant and product geometries were in fact local minima and not connected to the IRC path. For example, the geometry optimization was restarted using the final IRC geometries by removing the H and H 2 molecules. In all cases, the minima were found to be the same as obtained independently (as shown in Figures 2, 4). This suggests that the TS and the MEP are common to the ER and LH mechanisms. There are several possibilities to explain these results. First, the presence of an H atom or an H 2 molecule even at rather large distances has a strong enough influence to re-orient the H 2 CO and HCO molecules with respect to the surface. Nevertheless, the PES is rather flat too. Finally, the fact that the surrounding water molecules are frozen may have also some minor impact on the geometries. Table 1 shows the comparison between the ONIOM energies of the stationary points using the two levels of calculations [M05-2X/aug-cc-pVTZ and CBS/DLPNO-CCSD(T)]. The same trends found in gas phase occur, namely that the exothermicity decreases when going from M05-2X to CBS (−11.47 to12.79 kcal/mol for ΔE). On the other hand, the barrier height is reduced, from 8.28 to 6.75 kcal/mol. It is also possible to compare the values of these properties to those obtained in gas phase with the same methods. As can be seen in Table 1, the exothermicity is larger in gas phase for both methods, while the barrier height is only slightly larger. This would certainly result in gas phase rate constants being somehow larger at small temperatures than those obtained from the grain model. Thus, the ice model yields a noticeable stabilization of the gas-phase TS with respect to the effect on the minima. Song and Kästner (2017), at the PWB6K-D3/def2-TZVP level, found the same behavior going from gas phase to pure ASW surface, i.e., an augmentation of the barrier height by about 1 kcal/mol and a reduction of the exothermicity ( Table 1). This latter reduction is somewhat larger in Song and Kästner (2.41 kcal/mol) compared to the present result (1.01 kcal/mol). Another interesting quantity which can be obtained from the present calculations is the binding energy (BE) of the reactants and products to the ice surface [see also (Duflot et al., 2021)]. We define the BE of a molecule M on the icy grain as the positive quantity: where E (ice + M) is the energy of the adsorbed molecule M on the icy grain, E(M) is the energy of M in the gas phase and E(ice) is the energy of the bare icy grain. The calculations are performed with the QM/MM procedure detailed above, and include ZPE from DFT. For H 2 CO, BE is found to be 9.18 kcal/ mol at DFT level and 8.38 kcal/mol at CBS level. At the PWB6K-D3/def2-TZVP/MM level (without ZPE), Song and Kästner (2017) found a broad distribution of BE's on 81 structures obtained on pure ASW, between about 2.0 and 18.6 kcal/mol. The distribution was centered on BE's between about 8.0 and 9.9 kcal/mol. These authors selected one ASW structure upon which they calculated the rates with a BE of about 9.14 kcal/mol. Finally, our QM/MM BE value lies within accordance with their BE distribution. For HCO, the DFT and CBS values we computed are 8.08 and 7.37 kcal/mol, respectively. For this latter species, the CBS values can be compared to the largest BE's obtained on pure crystalline (7.08 kcal/mol) and amorphous (5.67 kcal/mol) ices (Duflot et al., 2021). However, it should be kept in mind that in the present work, the ice surface was not allowed to relax in reactivity studies, and that in ref (Duflot et al., 2021) the ice  Frontiers in Astronomy and Space Sciences | www.frontiersin.org March 2022 | Volume 9 | Article 807649 model is either pure ASW, or perfect crystalline Ih ice. Relaxation of the surface during the reaction mechanism could alter the BE even in the case of physisorption. in (Siaï et al., 2016)], due to the single point CCSD(T)//MP2 level, the present DFT and CBS V MEP profiles are more regular. The CBS V G a curve displays two small bumps at about s = 0 and s = 0.27 amu 1/2 bohr, separated by a small well at s = 0.15 amu 1/2 bohr. The inclusion of ZPE corrections does not change the position of the maximum of the curve which remains very close to s = 0 amu 1/2 bohr. The ΔZPE curve shows the evolution of the vibrational frequencies, using the TS as a reference: around the TS, i.e., between s ≈ 0.2 and s ≈ −0.3 amu 1/2 bohr, the variation remains close to 0. Then a steep variation occurs and rapidly leads to the asymptotic values corresponding to the formation of the formaldehyde CH bond (left side) and the formation of the H 2 bond (right side).

Kinetic Study: H 2 Formation Mechanisms
Eley-Rideal Mechanism Figure 7A displays the calculated rates for the Eley-Rideal mechanism, compared to the analytical functions fitted using Eq. 11 and the results from Song and Kästner (2017). It should be kept in mind that Song and Kästner used instanton theory (Zaverkin and Kästner, 2020) between 59 K and the so-called crossover temperature (Affleck, 1981) T c = 366 K and variational TST for higher temperatures up to 1000 K. Also, in Eq. 11, the β parameter was set to 1.0 in ref. (Song and Kästner, 2017). In the present work, β was optimized and the temperature spans a broader interval (15-2000 K). The fitted parameters for CVT rates are shown in Table 5 and compared to Song and Kästner. From Figure 7A, it is clear that Eq. 11 fits very well both ZCT and SCT rates. At very high temperatures, tunnelling effects are small (e.g., κ ≈ 1.1 at 2000 K for both ZCT and SCT). However, at room temperature, tunneling effects are still important with κ(ZCT) = 16.8 and κ(SCT) = 34.0. As recently pointed out by Richardson (2018), when the TS imaginary frequency is larger than 1300 cm −1 , deep-tunnelling still occurs at this temperature. This is indeed the case in the present work since the TS imaginary frequency is 1737 cm −1 ( Table 3). This gives a cross-over temperature of T c ≈ 400 K, which is very similar to that of Song and Kästner (2017). At low temperatures, the CVT/ ZCT rates appear to be several orders of magnitude lower than the CVT/SCT, which is a common trend (Truhlar and Kuppermann, 1971;Liu et al., 1993), due to the corner cutting effects (Allison and Truhlar, 1998). Our CVT/SCT curve is thus closer to Song and Kästner's one than the CVT/ZCT one. There remains however a difference of several orders of magnitude as can be seen in Table 6, between the present rates and the values provided by Song and Kästner for selected temperature values. The differences are due to the different barrier heights, IRC shape and theoretical models used to compute the rates. These differences are also reflected in Table 5: there is a large difference between the ZCT and SCT fitted parameters. Furhtermore, β is significantly different from 1 in our case.

Langmuir-Hinshelwood Mechanism
In Figure 7B, the calculated rate constants obtained for the LH mechanism are shown together with Song and Kästner (2017)'s Frontiers in Astronomy and Space Sciences | www.frontiersin.org March 2022 | Volume 9 | Article 807649 fitted rate function. As in the ER case, the CVT/ZCT and CVT/ SCT curves were fitted using Eq. 11. The optimized parameters are given in Table 5. Most of the trends observed in the ER case, such as the importance of tunnelling at high and low temperatures and the difference between the ZCT and SCT results, are also observed. The main difference comes from the comparison between the present CVT/SCT rate and Song and Kästner's one: both curves appear to be very close, especially at low temperature. This is reflected in the fitted parameters given in Table 5: although our α value is about 40 times larger, the other parameters remain very similar. In contrast to our ER results, β is now very close to 1. Table 6 also shows that the rates differ by a factor of around 3 at all selected temperatures. The similarity with Song and Kästner's results may be rationalized by the closeness of our LH barrier ( Table 4) compared to ref (Song and Kästner, 2017). Indeed, these authors calculated the LH rates using a complex with a barrier of 3030 K or 6.02 kcal/mol, virtually identical to the ER barrier. As shown in Table 4, our LH barrier is very close to this value (6.52 kcal/mol at the CBS level). Given the large rate coefficients predicted in the LH mechanism even at T = 15 K, it is anticipated that besides the exothermicity of the reaction process, the reaction mechanism is in favor of the formation of gaseous H 2 molecules. However, even though our ice model is intended to mimic ISM conditions, we would expect that diffusion of the incoming H (gas) atoms, so that they meet one CH 2 O on the grain, will be the rate limiting step.

CONCLUSION
In the present work, we have employed Density Functional Theory for geometry optimization, ZPE calculations, and CBS/ DLPNO-CCSD(T)//ONIOM/aug-cc-pVTZ single point energies methods on a large-scale system derived from a classical molecular dynamics trajectory. These calculations serve indeed as a feasibility assessment of previously benchmarked DFT (and to some extent, CVT) methods for the study of chemical reactions.
We study the abstraction reaction between an incoming gas phase hydrogen atom and another hydrogen atom belonging to a formaldehyde molecule at the surface of a mixed formaldehyde/water icy grain. The size of the system required the use of hybrid QM/MM techniques to obtain the detailed shape of the full reaction pathway including the effect of the environment. The entire minimum energy path was computed with CBS and M05-2X zero-point energy inclusion for each stationary point.
These results were then used to obtain the rate constants of this reaction between 15 and 2000 K for both ER and LH mechanisms using methods more sophisticated than the usual transition state theory, including contributions from variational and tunnelling effects. This work provides new results with respect to the choice of the model grain, i.e., doped formaldehyde-H 2 O amorphous ice, to the application of CBS model chemistry as opposed to pure DFT treatments, and to the use of VTST over a larger temperature range than what was reported previously. Our calculations show that both reaction models (ER and LH) exhibit exothermic pathways.
The present work focused on the H abstraction of H 2 CO and formation of H 2 , but the H 2 CO + H reaction has actually several possible products, such as H 3 CO, H 2 COH. Using the proposed methodology to characterize the latter products formation through reactions at the grain surface would definitively be of interest and complementary to the present work.
In addition, the reaction was fully characterized for one configuration that was arbitrarily chosen along a classical molecular dynamics trajectory. It can be expected that various binding sites may exist on the surface characterized by different environments around the surface CH 2 O molecule. To confirm the trend, a sampling should be carried out over a variety of binding sites. However, a full characterization of the reaction path, as presented here, has to be carried out for each binding environment, such an investigation certainly requires a substantial computational effort.
Finally, the reaction was characterized on a mixed grain, the reacting CH 2 O molecule forming part of the doped ice model. This work proposes a methodology to account for the presence of impurities in interstellar ices. Further efforts are needed to quantify more precisely the influence of these dopants on reaction rate constants both in terms of molecular or atomic diversity and concentration of the species.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
DD did some of the quantum chemical calculations, data curating and wrote parts of the manuscript. BK performed all DFT calculations, and reaction rates, and provided a first draft of the paper. MB prepared the QM/MM input structures for Gaussian calculations. MH performed the MD calculations. CT supervised the MD calculations and participated in the data curating.

ACKNOWLEDGMENTS
The authors thank the PCMI (Programme National de Physique et Chimie du Milieu Interstellaire) for support and the Centre de Ressources Informatiques (CRI) of the Université of Lille for providing computing time. MH acknowledges Dr C. Fotsing-Kwetche and Dr E. Michoulier for their help on the Molecular Dynamics part. MH acknowledges the Ministry of Higher Education and Scientific Research Tunisia for 2 exchange scholarships with Laboratoire de Physique des Lasers, Atomes et Molécules, UMR CNRS 8523, Université de Lille. DD and CT acknowledge the support of the OVERSEE and CAPPA grants, managed by the Agence Nationale de la Recherche under the frame programs Investissements d'Avenir ANR-10-LABX-005 and I-SITE ULNE/ANR-16-IDEX-0004 ULNE, respectively. They also thank the Région Hauts de France, the Ministre de l'Enseignement Supérieur et de la Recherche (CPER Climibio) and the European Fund for Regional Economic Development for their support. This work used HPC resources from GENCI-TGCC (Grant No. 2020-A0050801859).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fspas.2022.807649/ full#supplementary-material Supplementary Figure S1 | Radial distribution functions (A) between the water oxygens in case of doped ice and pure ice at T=110K; (B) between the carbon atoms of the H 2 CO molecules Supplementary Table S1 | Complete calculated reaction ER rate constants in cm 3 molecule -1 s -1 for H 2 formation on the grain model at CBS level.