Paper_Uracil_collision/paper_UAR.tex

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This is a (brief) model paper using the achemso class
%% The document class accepts keyval options, which should include
%% the target journal and optionally the manuscript type. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentclass[journal = jpcafh ,manuscript=article]{achemso}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Place any additional packages needed here.  Only include packages
%% which are essential, to avoid problems later. Do NOT use any
%% packages which require e-TeX (for example etoolbox): the e-TeX
%% extensions are not currently available on the ACS conversion
%% servers.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage[version=3]{mhchem} % Formula subscripts using \ce{}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% If issues arise when submitting your manuscript, you may want to
%% un-comment the next line.  This provides information on the
%% version of every file you have used.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%\listfiles

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Place any additional macros here.  Please use \newcommand* where
%% possible, and avoid layout-changing macros (which are not used
%% when typesetting).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand*\mycommand[1]{\texttt{\emph{#1}}}

%% Added by author
\usepackage{graphicx}
\usepackage{tikz}
\usepackage{xcolor}
\usepackage{pgfplots}
\usepackage{epstopdf}
\usepackage{multirow}
\usepackage{array}

\newcommand{\red}[1]{{\color{red}#1}}
\newcommand{\blue}[1]{{\color{blue}#1}}

\usetikzlibrary{calc,trees,positioning,arrows,chains,shapes.geometric,%
                decorations.pathreplacing,decorations.pathmorphing,shapes,%
                matrix,shapes.symbols,shadows}

\newcolumntype{C}[1]{>{\centering\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}

\graphicspath{{images/}}

\definecolor{cream}{RGB}{222,217,201}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Meta-data block
%% ---------------
%% Each author should be given as a separate \author command.
%%
%% Corresponding authors should have an e-mail given after the author
%% name as an \email command. Phone and fax numbers can be given
%% using \phone and \fax, respectively; this information is optional.
%%
%% The affiliation of authors is given after the authors; each
%% \affiliation command applies to all preceding authors not already
%% assigned an affiliation.
%%
%% The affiliation takes an option argument for the short name.  This
%% will typically be something like "University of Somewhere".
%%
%% The \altaffiliation macro should be used for new address, etc.
%% On the other hand, \alsoaffiliation is used on a per author basis
%% when authors are associated with multiple institutions.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\author{Kseniia A. Korchagina}
\affiliation[Universit\'{e} de Toulouse (UPS) and CNRS]{Laboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, Université\'e de Toulouse (UPS) and CNRS, 118 Route de Narbonne, F-31062 Toulouse, France}
\author{Fernand Spiegelman}
\affiliation[Universit\'{e} de Toulouse (UPS) and CNRS]{Laboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, Université\'e de Toulouse (UPS) and CNRS, 118 Route de Narbonne, F-31062 Toulouse, France}
\author{J\'er\^ome Cuny}
\affiliation[Universit\'{e} de Toulouse (UPS) and CNRS]{Laboratoire de Chimie et Physique Quantiques LCPQ/IRSAMC, Université\'e de Toulouse (UPS) and CNRS, 118 Route de Narbonne, F-31062 Toulouse, France}
\email{jerome.cuny@irsamc.ups-tlse.fr}
\phone{+33 (0)5 61 55 68 36}
\fax{+33 (0)5 61 55 60 65}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The document title should be given as usual. Some journals require
%% a running title from the author: this should be supplied as an
%% optional argument to \title.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title[An \textsf{achemso} demo]
{Molecular Dynamics Study of the Collision-Induced Reaction of H with CO on Small Water Clusters\footnote{Molecular Dynamics Study of the Collision-Induced Reaction of H with CO on Small Water Clusters}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Some journals require a list of abbreviations or keywords to be
%% supplied. These should be set up here, and will be printed after
%% the title and author information, if needed.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\abbreviations{PES,ISM,COSAC,FT-IR,MRCI+Q,CCD,CCSD(T),RCCSD(T),UCCSD(T),MP2, MBPT, BSSE, ZPE, DFT, SCC-DFTB, CM3, MD, MDPT}
\keywords{Molecular Dynamics, SCC-DFTB, HCO Radical, Collision-Induced Reaction, Water Clusters}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The manuscript does not need to include \maketitle, which is
%% executed automatically.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The "tocentry" environment can be used to create an entry for the
%% graphical table of contents. It is given here as some journals
%% require that it is printed as part of the abstract page. It will
%% be automatically moved as appropriate.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{tocentry}
\includegraphics[scale=0.1]{TOC.png}
\end{tocentry}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The abstract environment will automatically gobble the contents
%% if an abstract is not used by the target journal.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
The  successive hydrogenation of CO  is supposed to be the main mechanism leading 
to the formation of complex oxygenated species in the interstellar medium, possibly 
mediated by ice layers or ice grains. In order to simulate the dynamical influence 
of a water environment on the first step of the hydrogenation process, we achieve 
molecular dynamics simulations of the reactive collision of H with CO adsorbed on 
water clusters in the framework of the self-consistent-charge density functional 
based tight-binding approach (SCC-DFTB) to calculate Potential Energy Surfaces. The reaction 
probabilities and the reactive cross sections are determined for water cluster sizes 
up to ten water molecules. The collision results are analyzed in terms of different reaction pathways:
reactive  or non-reactive, sticking or desorption of the products or reactants. We show that the 
HCO radical, although potentially formed as an intermediate whatever the size of the 
water cluster, is significantly stabilized for cluster sizes larger than one water molecule 
and may remain adsorbed on water clusters with more than three molecules. This behavior 
is shown to be linked to the dissipation of the collision energy into vibrational 
excitation of the water cluster.
\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Start the main part of the manuscript here.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
In recent decades, the chemical composition of the Universe has been continuously 
investigated, essentially via absorption spectroscopy with background stars as lamps. 
It was established that all chemical elements and molecules in the InterStellar Medium 
(ISM) are concentrated in three main kind of clouds: diffuse, translucent and dense clouds. 
The chemical composition of these clouds is dominated by the abundance of hydrogen, 
the concentration of helium is about 10\%, and other elements such as carbon, nitrogen, 
and oxygen are also present in the ratio $\sim$10$^{-3}$-10$^{-4}$ of the hydrogen 
density \cite{Festou}. In addition to these elemental abundances, a large number of 
simple (H$_{2}$, CH, CH$^{+}$, CN, C$_{2}$, OH, CO, HCO, HCO$^{+}$, HCN, C$_{3}$) 
and more complex organic molecules (CH$_{3}$OH, C$_{2}$H$_{5}$OH, C$_{2}$H$_{5}$CN, 
CH$_{3}$COCH$_{3}$, CH$_{4}$, NH$_{3}$, H$_{2}$O, \textit{etc.}) were 
detected \cite{churchwell1,thompson1,schutte1,grim1,cordiner1,kaiser1,herbst1,herbst2}. 
Although most of those detected compounds can form through gas-phase processes,
some of them, in particular H$_{2}$, H$_{2}$O and CH$_{3}$OH are assumed to form
through chemical reactions between atoms and molecules at the surface of grains.\cite{Horn2004,Oberg2016}
The composition and structure of the outer layers of these grains can be of different nature: organic, 
ice-mineral mixture, pure water-ice and ice-dust mixture grains depending, for instance, 
on the incidence of UV radiation and cosmic rays\cite{schulz1}, or the atomic densities of
the environment where they are present.\cite{Taquet2012} In 2015, a unique set of experimental data on the chemical 
composition of the surface of the 67P/Churyumov-Gerasimenko comet were obtained through 
various instruments of the Philae module \cite{wright1,goesmann1,spohn1}. For instance, 
mass spectroscopy measurements performed by the Ptolemy instrument detected H$_{2}$O, 
CO$_{2}$ and CO as main volatile species in a ratio of 10:2:$>$1 \cite{wright1}. 
The COmetary SAmpling and Composition (COSAC) experiment also provided a picture of 
the organic composition of the comet surface.\cite{goesmann1} Sixteen molecules, 
belonging to two main molecular groups, were identified by the COSAC instrument. 
The first group includes the H$_{2}$O and CO molecules and the subsequent oxygen-containing 
organic species such as alcohols and carbonyls. The second group includes CH$_{4}$ 
and NH$_{3}$ (although not confidently identified) that lead to the formation of 
nitrogen-containing organic species (amines, nitriles, amides, and isocyanates) \cite{goesmann1}. 
These results provide a rather accurate picture of the 67P chemical composition 
(and potentially of other comets) but no accurate information is provided on the 
mechanisms that lead to the formation of these complex organic compounds in the 
ISM \cite{loomis1,kalvans1}. This is an important question to address for both 
theoreticians and experimentalists in order to strengthen our understanding of 
the ISM chemistry.

Currently, it is supposed that the main mechanism leading to the oxygenated 
compound formation in grain mantles and comets is a process of successive 
hydrogenation of CO, occurring by the consecutive addition of hydrogen 
atoms:\cite{book1,tielens1,crovisier1,watanabe2,watanabe3,hidaka1}

\begin{center}
CO $\xrightarrow{\text{H}}$ HCO $\xrightarrow{\text{H}}$ H$_{2}$CO $\xrightarrow{\text{H}}$ H$_{3}$CO $\xrightarrow{\text{H}}$ CH$_{3}$OH
\end{center}

This mechanism has been the subject of various experimental studies. For instance, Hiraoka
and co-workers studied the reaction of H atoms with a solid CO thin film in the 10-25~K
temperature range.\cite{hiraoka1,hiraoka2}  From these experiments, the authors
concluded that -(i)- H atoms do not diffuse into the CO matrix, the reaction proceeds only
at the surface -(ii)- the primary reaction product H$_{2}$CO is prone to polymerization -(iii)- the
rate constant for the first and second steps of this reaction is small at cryogenic temperature -(iv)-
formaldehyde and methanol were formed. Later, Pirim \textit{et al.}  studied the same reaction
in two different ways, by simple hydrogenation of a CO surface and by co-injection of CO molecules
and H atoms \cite{pirim1}. In the former case, nothing was formed at 3~K whereas H$_{2}$CO and
CH$_{3}$OH were detected at 10~K by Fourier Transform InfraRed (FT-IR) spectroscopy. In contrast,
the co-injection at 10~K leads to the formation of both the HCO and H$_{3}$CO radicals as major
products. However, as mentioned above, the   hydrogenation of CO in the ISM is likely to
occur at the surface of ice-coated grains. To understand such a situation, some groups also studied
the H addition on CO in a H$_{2}$O-CO ice \cite{watanabe1,watanabe2,pirim2}. They showed that
water molecules play two important roles in this process. On the one hand, they exhibit a catalytic
role by helping to overcome the activation barriers. On the other hand, water molecules create new
chemical pathways that enhance the reactivity.

To complement these experimental measurements, theoretical studies were conducted on the different 
steps of the sequential hydrogenation reaction of carbon monoxide, isolated or in the presence of 
water molecules \cite{woon1,cao1,woon2,rimola1,peters1,Peters2013a}. In particular, the geometric and 
energetic characteristics of reactants, transition states and products involved in this mechanism 
have been the subject of various calculations at high levels of theory. For instance, Woon observed 
that, at the QCISD$^{*}$ level of theory, the barrier of the H + H$_{2}$CO $\longrightarrow$ 
CH$_{3}$O reaction is reduced from 5.04~kcal.mol$^{-1}$ for isolated CO to 4.18~kcal.mol$^{-1}$ 
when adding four water molecules \cite{woon1}. The intermediates of this multi-step mechanism, 
such as HCO and H$_{2}$CO, may also undergo isomerization under certain conditions and were 
thus also investigated \cite{zanchet1,schreiner1,peters1}.

However, it is certainly the first step of the sequential hydrogenation of CO that has deserved 
the largest attention in the literature as it is a crucial step towards the formation of complex 
oxygenated species \cite{woon1,cao1,woon2,rimola1,peters1,Peters2013a,adams1,bitter1,werner1,werner2,cho1,romanowski}.
For instance, Woon compared the structures and frequencies of CO and HCO as well as the energetics 
of the H + CO $\longrightarrow$ HCO reaction (without water) with various high-accuracy schemes 
such as Restricted Coupled Cluster RCCSD(T) or Davidson-corrected internally-contracted
MRCI+Q methods \cite{woon2}.  The author showed that at the aug-cc-pVTZ level, MRCI+Q and
RCCSD(T) methods with ZPE correction  slightly overestimate barrier heights, 4.1 and 4.2~kcal.mol$^{-1}$
respectively, in comparison  with the experimental value 2.0$\pm$0.4 kcal.mol$^{-1}$. In contrast,
the MRCI+Q and RCCSD(T) results for the reaction ergicity, -13.7 and -13.8 kcal.mol$^{-1}$, respectively,
are in a very good agreement with the experimental value of -14 kcal.mol$^{-1}$. Those calculations
were recently revisited by Peters \textit{et al.} using state-of-the-art multireference \textit{ab initio}
methods to study the energetics of the HCO/DCO formation and dissociation processes.\cite{Peters2013a}
In addition to their theoretical interest, the results of those and other \textit{ab initio} calculations 
can be further used as input values for larger scale numerical simulations. For instance, the GRAINOBLE model
was used by Rimola  \textit{et al.} \cite{rimola1} to describe the distribution of the H$_{2}$CO and CH$_{3}$OH 
ice abundances. It leads to abundance values that are in a good agreement  with the experimental data.
Although extremely informative, high-precision \textit{ab initio}  wavefunction type calculations are only
achievable for very small systems. In a similar way, Density Functional  Theory (DFT) can describe larger species
and has been applied to systems containing up to 32 water molecules  but only in a static framework.\cite{rimola1} 
Although a good description of the Potential Energy Surface (PES) of the various species 
involved in the mechanism is of primary importance, a dynamical simulation  of the hydrogenation 
process at the molecular level is also desirable. Indeed, size and temperature effects of the water
substrate, influence of its morphology  and diffusion of species at the surface of nanograins can
be important  contributing factors that can hardly be included in highly accurate quantum chemical 
calculations although they can have a strong impact. Furthermore, although reactions can occur at the 
surface of grains between diffusing species \cite{Oberg2016}, they are also likely to occur during 
collisions between molecules. This can hardly be described by static quantum chemical calculations 
only and thus a more complete understanding of such a mechanism require the use of Molecular 
Dynamics (MD) simulations. To the best of our knowledge, no MD study has been conducted to 
explicitly simulate the reaction of H with CO in the presence of water molecules.

In the present contribution, we present a MD study of the collision of H with CO at the surface 
of water clusters, \textit{i.e.} CO-(H$_{2}$O)$_{n}$ (n=0-10), aiming at a statistical sampling of
the various possible pathways likely to occur between the two species. Such simulations are
made possible by the description of the PES at the Self-Consistent-Charge Density Functional 
based Tight-Binding (SCC-DFTB) level of theory that allows to describe bond-forming and
bond-breaking at a limited computational cost.  The outline of the article is as follows: the
computational model is described in Section II,  the results of the simulations are presented in
Section III and the main outcomes and perspectives are summarized in the conclusion.

\section{Computational Methods} \label{Comput_meth}

\textbf{Computation of the PES.}
We use the SCC-DFTB \cite{elstner,koskinen,frauenheim} approach implemented in the deMonNano code \cite{heine1} 
to describe the PES of the CO-(H$_{2}$O)$_{n}$ clusters. Within the SCC-DFTB formalism, the electronic
energy is given by the following equation:
%
\begin{eqnarray}\label{enr}
E^{SCC-DFTB} = \sum\limits^{occ} \langle \psi_i | \hat{H_0} | \psi_i \rangle + \sum\limits_{\alpha \beta} U_{\alpha \beta} (R_{\alpha \beta})
+ \frac{1}{2} \sum\limits_{\alpha \beta} \Delta q_{\alpha} \Delta q_{\beta} \gamma_{\alpha \beta} 
- \sum\limits_{\alpha \beta} f_{damp} \frac{C_{6}^{\alpha \beta}}{R_{\alpha \beta}^{6}}
\end{eqnarray}
%
where the 1$^{\text{st}}$ term is a tight-binding term defined from parametrized integrals and 
the 2$^{\text{nd}}$ term is a repulsive interaction expressed as a sum over all atomic pairs. 
In the present study, we used the mio-set for Slater-Koster integrals \cite{elstner}. 
The 3$^{\text{rd}}$ term is the second-order term of the Taylor expansion expressed as a function 
of the atomic charge fluctuations $\Delta q_\alpha$ and the 4$^{\text{th}}$ term describes the 
London dispersion interaction. Rapacioli and co-workers proposed to improve the description of 
the electrostatic interaction in molecular systems by replacing the original Mulliken charges by 
the Class IV - Charge Model 3 (CM3) charges \cite{DFTB_CM3,rapacioli2}, defined as:
%
\begin{equation}\label{cm3}
q_{\alpha}^{CM_3} = q_{\alpha}^{Mull} + \sum \limits_{\alpha' \neq \alpha}^{atoms} [D_{Z_{\alpha}Z_{\alpha'}}B_{\kappa \kappa'} + C_{Z_{\alpha}Z_{\alpha'}}B^{2}_{\kappa \kappa'}]
\end{equation}
%
where $B_{\kappa \kappa'}$ is the Mayer's bond order whereas $C_{Z_{\alpha}Z_{\alpha'}}$ and
$D_{Z_{\alpha}Z_{\alpha'}}$ are empirical parameters to define. We used the D$_{\text{OH}}$=0.129
value previously proposed by Simon and Spiegelman in their study  on water clusters \cite{simon1,simon2,simon3}. 
We fitted the D$_{\text{CO}}$ parameter to reproduce the experimental dipole moment of CO which
leads to D$_{\text{OC}}$=0.012. We set D$_{\text{CH}}$=0.0.

A difficulty in the description of the H+CO reaction lies in the fact that it involves, 
at small distances, a potential crossing between the  $^{2}\Pi$ electronic ground-state 
(which dissociates to ground state products H($^{2}S$) + CO($^{1}\Sigma^{+}$)) and a 
$^{2}\Sigma^{+}$ excited state of the formyl radical correlated with H($^2S$)+CO($^3\Pi$).
Therefore, at the intersection of the electronic states, the treatment of the problem with 
a mean-field single determinant scheme arises, and self-consistency convergence problems 
occur. To solve this problem, we used a Fermi-Dirac orbital occupation defined by an
electronic temperature of 1000~K which allows to achieve a continuous switch from
one state to the other in the  near vicinity of the crossing.

\textbf{Exploration of the PES.}
Before performing the collisional trajectories between H and CO, we first optimized the 
geometry of the considered CO-(H$_{2}$O)$_{n}$ clusters in order to start our trajectories 
as close as possible to the equilibrium configurations at low temperature. To explore the 
PES of the clusters in an exhaustive way, we used the Molecular  Dynamics Parallel-Tempering
(MDPT) algorithm \cite{sugita1,sugita2,earl1}, which allows for replica 
exchanges between trajectories at different temperatures. This scheme increases 
the ergodicity of the MD simulations and thus speeds up the exploration of the PES at a determined 
temperature. We used a temperature range going from 20 to 320~K by steps of 5~K which correspond 
to 60 distinct temperatures. All the MD trajectories were 4~ns long with a timestep of 0.2 fs. 
The PT replica exchanges were attempted every 400 fs. In order to achieve canonical simulations, 
we used a Nos\'e-Hoover chain of five thermostats defined by a unique  frequency of 800 cm$^{-1}$ \cite{nose,hoover}.

In order to find the lowest-energy configurations on each PES, local geometry optimizations 
of the CO-(H$_{2}$O)$_{n}$ clusters were subsequently performed using the following procedure: 
for one out of  four temperatures, \textit{i.e.} every 20~K, one thousand different geometries 
were periodically selected and locally optimized using a conjugate gradient algorithm. This lead 
to a total of 15000 geometry relaxations per CO-(H$_{2}$O)$_{n}$ cluster,  from which  the
lowest energy  one was retained.

\textbf{MD Collision Trajectories.}
The first step of the MD simulation consisted in the sampling of the initial conditions, namely 
-(i)- the initial positions and velocities of the CO-(H$_{2}$O)$_{n}$ clusters, -(ii)- the angular 
orientations of the clusters and -(iii)- the impact parameter of the collision. In order to sample
the initial positions and velocities, we achieved thermalisation of the CO-(H$_{2}$O)$_{n}$ 
clusters at 70~K via a 200~ps long MD simulation in the canonical ensemble using the previously
obtained lowest-energy 
configuration as initial geometry. The last geometries and corresponding velocities were then taken 
as initial conditions for subsequent collisional  trajectories. These initial positions of the
CO-(H$_{2}$O)$_{n}$ clusters were further evenly rotated along the three Cartesian axes to obtain
64 initial angular conditions. For each orientation, the impact parameters 
defining the initial position of the colliding hydrogen atom were randomly generated in a disk 
of radius R (R$>$R$_{\text{cluster}}$) centered at 10 \AA \ from the center of mass of
the cluster. For the CO-(H$_{2}$O)$_{n}$ clusters with n=0-5, 53 different impact parameters
were generated per angular orientation leading to a sampling of 3392 initial conditions.
In the case of CO-(H$_{2}$O)$_{10}$, 2000 initial positions of the hydrogen atom were generated
opposite to the surface containing the CO molecule in order to describe quasi-frontal collisions only.
The initial velocity of the hydrogen atom  was set to 0.01 \AA.fs$^{-1}$ which is consistent with
the temperature of the cluster at 70~K.  The time length of each collision trajectory was set equal to 10~ps.

\textbf{Wavefunction and DFT Calculations.}
In order to establish reference equilibrium geometries and intermolecular interaction energies 
of the CO-H$_{2}$O isomers, we performed \textit{ab initio} MP2 and CCSD(T) calculations in 
combination with the Pople-style 6-311++G(d,p) basis set \cite{krishnan1} and the aug-cc-pVTZ 
basis set of Dunning and co-workers \cite{dunning1,kendall1}, respectively. To further check 
the performances of the SCC-DFTB approach, we also carried out DFT calculations using the 
B3LYP,\cite{becke1,lee1,vosko1} B3LYP-D3\cite{Grimme2010} and B97D\cite{Grimme2006}
exchange-correlation functional in combination with the 6-311++G(d,p) basis set. The latter
functional was tested as it was shown by Peters \textit{et al.} to provide a satisfactory value
for the formation barrier of the H + CO $\rightleftharpoons$ HCO reaction in vacuum.\cite{Peters2013a}
Basis set superposition  errors (BSSE) were taken into account using the counterpoise method of Boys
and Bernardi.\cite{Boys2002} All DFT and wavefunction calculations were performed with the Gaussian
09 package \cite{g09}.
 
All the binding energies between CO and the water clusters discussed in the text were defined as the energy of the
relaxed CO-(H$_{2}$O)$_{n}$ complex minus the energy of the water cluster minus the energy of CO both taken in
their geometry in the optimized complex. In the same way, formation energies for the H+CO-(H$_{2}$O)$_{n}$ reaction
were defined as the energy of the relaxed HCO-(H$_{2}$O)$_{n}$ complex minus the energy of one hydrogen minus the
energy CO-(H$_{2}$O)$_{n}$ considering its geometry in the optimized complex.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and Discussion} \label{resul_disc}

\subsection{Validation of SCC-DFTB Potential.}

The SCC-DFTB method is a parametrized approach. The D$_{\text{OH}}$ parameter involved in the 
definition of the CM3 charges for the O--H pair was shown by Simon \textit{et al.} to yield 
good results for the description of water clusters.\cite{simon1,simon2,simon3}
For instance, in the water dimer, it provides an oxygen-oxygen distance of 2.92 \AA \ to be 
compared with the experimental value of 2.98 \AA \cite{odutola}. The validity of the 
D$_{\text{CO}}$ and D$_{\text{CH}}$ parameters determined in the present work had to be 
checked. For this purpose, we compared  the structural and energetic characteristics of 
the two well-known CO-H$_{2}$O isomers to fully \textit{ab initio} calculations. 
Figure~\ref{str} shows the structure of these two CO-H$_{2}$O isomers. The first one 
(C-structure) corresponds to a bonding  between an H atom of water and the C atom of 
CO while the second one (O-structure) corresponds to the interaction between an H atom 
of water and the O atom of CO.

\begin{figure}[h!]
\centering
\includegraphics[scale=0.40]{H2O_CO_OC.png}
\caption{MP2/6-311++G(d,p) geometries of the global (C-structure) and local  (O-structure)  minima of
CO-H$_{2}$O.} \label{str}
\end{figure}

\begin{table}[h]
\small
\centering
\caption{Equilibrium  geometries and intermolecular interaction energies of the C- and O-structures of
CO-H$_{2}$O obtained from SCC-DFTB, DFT/B3LYP, DFT/B3LYP-D3, DFT/B97D and MP2 with the 6-311++G(d,p)
basis set and CCSD(T) with the aug-cc-pVTZ basis set.}
\label{r_e}
\begin{tabular}{|c|c|c|c|c|}
\hline
\multirow{2}{*}{Method}  & \multicolumn{2}{c|}{C-structure} & \multicolumn{2}{c|}{O-structure} \\ \cline{2-5}
                         & $r_{min}$, \AA & \begin{tabular}[c]{@{}c@{}}E$_{\text{int.}}$,\\ kcal.mol$^{-1}$\end{tabular}    & $r_{min}$, \AA & \begin{tabular}[c]{@{}c@{}}E$_{\text{int.}}$,\\ kcal.mol$^{-1}$\end{tabular} \\ \hline
SCC-DFTB   & 2.20           & -0.72                     & 2.18           & -0.44      \\ \hline
MP2        & 2.45           & -1.84                     & 2.41           & -1.02      \\ \hline
MP2/BSSE   & 2.52           & -1.51                     & 2.53           & -0.61      \\ \hline
B3LYP      & 2.44           & -1.49                     & 2.41           & -0.79      \\ \hline
B3LYP/BSSE & 2.42           & -1.36                     & 2.40           & -0.66      \\ \hline
B3LYP-D3          & 2.38           & -2.03                     & 2.33           & -1.47      \\ \hline
B3LYP-D3/BSSE & 2.40           & -1.89                     & 2.36           & -1.28      \\ \hline
B97D           & 2.51           & -1.82                     & 2.50           & -1.06      \\ \hline
B97D/BSSE   & 2.53          & -1.68                     & 2.52           & -0.92      \\ \hline
CCSD(T)$^a$    & -              & -2.05                     & -              & -1.21      \\ \hline
CCSD(T)/BSSE$^b$ & -            & -1.65                     & -              & -0.90      \\ \hline
\end{tabular}\\
$^a$ CCSD(T) energy calculations were performed using the MP2 optimized geometry\\ 
$^b$ CCSD(T)/BSSE energy calculations were performed using the MP2/BSSE optimized geometry\\ 
\end{table}

The structural and energetic characteristics of these dimers are given in Table~\ref{r_e}. 
The interaction distances calculated by the SCC-DFTB method are 2.20 and 2.18 \AA~for r(H$_{2}$-C) 
and r(H$_{2}$-O$_{2}$), respectively, somewhat smaller  than the data obtained from MP2/6-311++G(d,p) 
calculations, namely 2.45 and 2.41 \AA, respectively. From an energy standpoint, the CO-H$_{2}$O is 
a very weakly bound complex (both C- and O-structures), where the interaction energy includes about 
40 \% of electrostatic, 35 \% of induction and 25 \% of dispersion contributions \cite{vilela1,wheatley1,yaron1,sadlej1}.
Referring to such a very small value of the interaction energy, an accurate description using 
density functional or semi-empirical methods is a very difficult task. Thus, although the SCC-DFTB 
correctly predicts the energetic ordering of the two isomers, the interaction energies are somewhat 
underestimated, even-though it  appears from  Table~\ref{r_e} that the exact value of the interaction 
energy is highly  dependent on the method, basis set, and applied corrections. Indeed, the three
DFT calculations (B3LYP, B3LYP-D3 and B97D) leads to quite different values for the interaction energy.
Interestingly, the B97D functional in combination with the 6-311++G(d,p) basis set provides geometries
that are very close to MP2/BSSE structures, and interaction energies almost equal to CCSD(T)/BSSE ones.
This confirms the good performances of B97D to describe carbon monoxide.\cite{Peters2013a}
From the Table~\ref{r_e}, the SCC-DFTB results for CO-H$_{2}$O are thus qualitatively similar to \textit{ab initio}
data and one can expect that the observed differences will not affect too seriously the present dynamical simulations.
As shown in the Table~\ref{e_bind}, with increasing cluster size, the SCC-DFTB binding energies converge rather
quickly from 0.72 kcal.mol$^{-1}$ to 1.00 kcal.mol$^{-1}$ and reach convergence at n=4-5. A similar behavior is
obtained at the B97D and MP2 levels of theory while B3LYP-D3 values display larger fluctuations. Consequently,
although the SCC-DFTB binding energies are still underestimated for larger clusters with respect to DFT and MP2 values,
the rapid convergence trend is expected to be reliable and certainly important for the reactive collisional behavior
discussed in the next section.

\subsection{Results}
\begin{figure}[h!]
\centering
\includegraphics[scale=0.45]{co_h2o_figures.png}
\caption{Geometries of the global SCC-DFTB minima of the
CO-(H$_{2}$O)$_{n}$ (n=1--5,10) clusters.} \label{str_2}
\end{figure}

\begin{table}[]
\centering
\caption{Binding energies between the CO (second column)  and HCO (third column) molecules and (H$_{2}$O)$_{n}$
clusters obtained with SCC-DFTB as well as B3LYP-D3, B97D and MP2 in combination with the 6-311++G(d,p)
basis set. \textit{Ab initio} values are reported including BSSE corrections.
SCC-DFTB HCO formation energy for H+CO-(H$_{2}$O)$_{n}$ along with B97D results
obtained by Peters \textit{et al.}\cite{} (Fourth column). All energies are in kcal.mol$^{-1}$.}
\label{e_bind}
\begin{tabular}{|c|cccc|cccc|c|}
\hline
  & \multicolumn{4}{c|}{CO-(H$_{2}$O)$_{n}$} & \multicolumn{4}{c|}{HCO-(H$_{2}$O)$_{n}$} &  H+CO-(H$_{2}$O)$_{n}$ \\ 
 \multirow{-2}{*}{n}   &  DFTB & B3LYP-D3 & B97D   & MP2    & DFTB    & B3LYP-D3 & B97D   & MP2     & DFTB \\ \hline
0                                &    -     &        -        &        -  &       -  &     -      &    -            &   -        &    -     &  -31.18 (\textit{-25.69}$^a$)       \\ \hline
1                                & -0.72 & -1.89        & -1.68  &-1.51  &-5.32    &-3.38         &-2.71    &-2.46  & - 26.25 \\ \hline
2                                & -0.93 & -3.62        & -2.86  &-2.52  & -10.86 &-8.22         &-6.63    &-5.50  & - 33.43 \\ \hline
3                                & -1.18  & -2.30       &-1.89   & -1.57 & -7.69   &-9.44         &-8.75    &-6.35  & - 27.95 (\textit{-26.76}$^a$) \\ \hline
4                                & -1.00  & -2.43       & -1.74  &-1.45  & -11.96 &-10.21       &-8.58    &-6.84  & - 26.22 \\ \hline
5                                & -1.00  & -1.99       & -1.78  & -1.53 &-4.77    &-5.87         &-4.62    &-3.47  & - 27.68 (\textit{-26.68}$^a$) \\ \hline
10                              & -1.01   &-1.88       &-1.92   & -1.81 &-3.54    & -3.93        &-2.84    &-1.26   & - 26.57 \\ \hline
\end{tabular} \\
$^a$ B97D results of Peters \textit{et al.}\cite{Peters2013a,Phillip2013}\\ 
\end{table}

To ensure that statistical convergence is reached, we checked the distributions of the initial
positions of the hydrogen atom with respect to the cluster in the 3392 starting configurations.
For visualisation, the distribution maps of the initial positions of hydrogen for CO-(H$_{2}$O)$_{n}$ (n=1,2,3)
are displayed in Figure~\ref{maps}. Similar data are obtained for CO-(H$_{2}$O)$_{n}$ (n=4,5).
From those pictures, it appears that the input geometries provide a reasonably uniform distribution
of hydrogen projectiles around the clusters. 

\begin{figure*}
\centering
\includegraphics[scale=0.45]{Maps_CO_H2O.png}
\caption{Distribution of the initial positions of the incident hydrogen atom around the CO-(H$_{2}$O)$_{n}$
         (n=1,2,3) clusters.} \label{maps}
\end{figure*}

For each cluster size, we investigated and analyzed all trajectories according to seven 
scenari (also called variants) characterizing the issue of the collision. A schematic 
description of these mechanisms is shown in Figure~\ref{Variants}.

\begin{figure*}
\centering
\includegraphics[scale=1.2]{Variants.png}
\caption{Possible pathways (variants) characterizing the issue of the H + CO collision.} \label{Variants}
\end{figure*}

Variant 1 corresponds to the case where  the two following conditions are fulfilled simultaneously:
-(i)- during the collision, the reaction occurred and the HCO radical was formed, 
and -(ii)- this radical remained connected to the cluster till the end of the simulation. 
Geometric characteristics of the formyl radical such as R$_{\text{CO}}$ and R$_{\text{CH}}$ 
were chosen to monitor the formation process. Recent studies by Adams and Purvis 
provide the accurate equilibrium bond lengths of HCO using Many-Body Perturbation Theory 
(MBPT) or Coupled-Cluster with Doubles excitations (CCD) calculations. The CCD values are 
R$_{\text{CO}}$=1.188 \AA \ and R$_{\text{CH}}$=1.111 \AA~ \cite{adams1,marenich1}. 
The present SCC-DFTB characteristics are R$_{\text{CO}}$=1.167 \AA \ and R$_{\text{CH}}$=1.197 \AA. 
The criterion for HCO formation during the dynamics was thus chosen as follows: if during the simulation 
the distances were fluctuating in the range  0.9-1.3 \AA~for C-O and  1.0-1.5 \AA~for C-H, 
we assumed that HCO was formed. To ensure that the radical remains bonded to the surface, 
we calculated the distances between the hydrogen atom of HCO and all oxygen atoms of the water 
cluster (R$_{\text{H--O}}$), and the distances between the carbon atom of HCO and all oxygen atoms 
of the water cluster (R$_{\text{C--O}}$). If at least one value of the R$_{\text{H--O}}$ 
distance and one value of the  R$_{\text{C--O}}$ distance were less than 4 \AA~at a given 
MD step, the radical was assumed to be associated with the water cluster.

The second variant  describes the situation where the HCO radical was formed, but then desorbed from the
cluster  (`Variant 2' in the Figure~\ref{Variants}). Similarly to CO, the binding energy between the formyl radical
and a water molecule is very small.Indeed, Cao \textit{et al.} calculated at the  UCCSD(T)/aug-cc-pVTZ
level of theory, corrected for BSSE and vibrational zero-point energies, the interaction energies of three stable
isomers of HCO--H$_{2}$O. They obtained the three following values: -0.83, -1.70 and -1.61~kcal.mol$^{-1}$.\cite{cao1}
The corresponding values without vibrational zero-point energies are -2.24, -2.62 and -3.35~kcal.mol$^{-1}$.
This weak interaction energy confirms that Variant 2 is likely to occur if the incident hydrogen atom has enough
kinetic energy. It should also be noted that, similarly to CO, the HCO--(H$_{2}$O)$_{n}$ binding energies
are rather low and quite sensitive to the computational details. This is illustrated in the Table~\ref{e_bind} for
different DFT and MP2 calculations. Although the binding energies for HCO--(H$_{2}$O) and HCO--(H$_{2}$O)$_{2}$
are somewhat overestimated, the DFTB values for larger species fall in the range of the DFT and MP2 values.
The next variant (`Variant 3' in the Figure~\ref{Variants})
corresponds to the  case where the reaction between H and CO did not occur while both reactants
remained stuck to the water cluster. The group of variants 4-6 may 
occur whenever the kinetic energy of the incident  hydrogen is above the energy required for 
the formation of the radical. They can be briefly described as follows: Variant 4: the HCO 
radical was not formed and H  desorbed from the cluster while  CO remained  adsorbed; 
Variant 5: the HCO radical was not formed, H was adsorbed while CO desorbed; Variant 6:
 the HCO radical was not formed, both H and CO  desorbed. It is worth pointing out that
those three variants can be encountered despite the transient formation of the HCO radical
if it dissociates during the trajectory. The last mechanism in our scheme  (`Variant 7' in the Figure~\ref{Variants})
is the situation where HCO is formed during some time of the simulation but finally dissociate
and both reactants remain bonded to the water cluster. The likelihood of each variant for clusters
CO-(H$_{2}$O)$_{n}$ (n=0-5,10) are listed in the Table~\ref{result_ratio}.

\begin{table*}
\small
\centering
\caption{Probability (in percent) of each collisional pathway as a function of the cluster size.}
\label{result_ratio}
\begin{tabular}[t]{C{2.1cm} C{1.4cm} C{1.4cm} C{1.4cm} C{1.4cm} C{1.4cm} C{1.4cm} C{1.4cm}}
\hline \hline
\small{System}  & \small{V.1 (\%)}  & \small{V.2 (\%)}  & \small{V.3 (\%)} & \small{V.4 (\%)} & \small{V.5 (\%)} & \small{V.6 (\%)} & \small{V.7 (\%)} \\
\hline
\small{CO}                    &  0.0   & 0.0   & 0.0   & 0.0  & 0.0   & 0.0   &  \textbf{35.3} \\
\small{CO-H$_{2}$O}           &  0.0  & 3.1   & 0.0  & \textbf{63.2}   & 9.1   & 10.2   & 7.5  \\
\small{CO-(H$_{2}$O)$_{2}$}   &  2.3   & \textbf{58.1}   & 0.3   & 5.0   & 9.7   & 0.0   & \textbf{24.6}  \\
\small{CO-(H$_{2}$O)$_{3}$}   &  \textbf{10.1}   & \textbf{50.9} & 1.3  & 0.9  & \textbf{31.4}   & 0.0   & 5.4 \\
\small{CO-(H$_{2}$O)$_{4}$}   &  \textbf{58.0}   & 0.0 & \textbf{35.8}  & 0.9  & \textbf{0.0}   & 0.0  & 5.4 \\
\small{CO-(H$_{2}$O)$_{5}$}   &  \textbf{61.4}   & 0.0 & 7.4  & 1.1  & \textbf{26.9}     & 0.1   & 0.8 \\
\small{CO-(H$_{2}$O)$_{10}$}   &  \textbf{34.8}  & 5.0 & \textbf{32.4}  & 3.0  & \textbf{20.8}     & 0.0   & 4.0 \\
\hline \hline
\end{tabular}
\end{table*}

In the gas phase  conditions, \textit{i.e} without water molecule, the formation of the HCO radical with a favorable
geometry occurred  in 35.3 \% of trajectories, while in other cases, the hydrogen atom flew away
from the collision  area. However, even in the favorable cases, the HCO radical always dissociated
at some point of the simulations. Indeed, after a successful reaction, the kinetic energy of the incident
hydrogen distributes into the few vibrational modes of the radical which leads to its dissociation.
To support this point and to highlight the role of the water cluster, Figure~\ref{r_ch} displays
a typical example of the time-evolution of the C--H distance for the HCO and HCO-(H$_{2}$O)$_{5}$
species. Immediately after the radical formation, the fluctuations of the C--H distance in the two
systems are of the same magnitude. However, after $\sim$1~ps, those fluctuations significantly
decrease for HCO-(H$_{2}$O)$_{5}$ highlighting the dissipation of kinetic energy towards the water
molecules. This leads to a stable HCO molecule. In contrast, the fluctuations are constant in the 
simulation of pure HCO which can result in its dissociation.

\begin{figure}[h!]
\centering
\includegraphics[scale=0.20]{rch_fluctuations.png}
\caption{Time-evolution of  the  C--H distance for HCO-(H$_{2}$O)$_{5}$ (top) and HCO (bottom)
immediately after the radical formation. $\Delta_{R_{C-H}}$ represents the amplitude of the
distance fluctuations.} \label{r_ch}
\end{figure}

The same behavior is observed for  CO-(H$_{2}$O). Indeed, $\sim$12 \% of the trajectories
lead to the formation of HCO although it is stable in only $\sim$3 \% of them. Variant
4, \textit{i.e.} with desorption of the hydrogen from the water molecule, is an important
pathway for this species as it encompasses $\sim$80 \% of the trajectories. For CO-(H$_{2}$O)$_{2}$
and CO-(H$_{2}$O)$_{3}$, the HCO radical is obtained in more than 50 \% of cases. Due to the
small binding energy of HCO with water molecules (see discussion above), most of the
successfully formed radicals  desorb from the cluster. This is demonstrated
by the predominance of the Variant 2 pathway. It is worth pointing out that the exact amount of
this latter pathway would be highly influenced by the level of theory used to describe the PES of
those systems. CO-(H$_{2}$O)$_{2}$ displays a significant amount
of dissociated radicals (24.6 \%) which supposes that two water molecules is not enough to
accommodate the excess kinetic energy of the hydrogen. In both CO-(H$_{2}$O)$_{2}$ and
CO-(H$_{2}$O)$_{3}$, only a few trajectories lead to a stable HCO radical adsorbed on the
water cluster, 2.3 and 10.1 \%, respectively.

Things are completely different for CO-(H$_{2}$O)$_{4}$ and CO-(H$_{2}$O)$_{5}$. Indeed,
there is no desorption of HCO from the cluster and only a few cases of HCO dissociation, 
5.4 and 0.8 \%, respectively. The majority of the simulations lead to a stable HCO radical
adsorbed on the water molecules. This allows us to assume that the second step of the
successive CO hydrogenation is likely to occur on complexes composed of four or more water
molecules. The main difference between those two species is the amount of Variant 3 and 5
pathways that seems to be opposed. In particular, CO-(H$_{2}$O)$_{4}$ displays a 0.0 \%
amount of Variant 5, in contrast to the other species.  This could likely be attributed to the
planar square conformation of (H$_{2}$O)$_{4}$.
Surprisingly, despite performing quasi-frontal collisions only, for CO-(H$_{2}$O)$_{10}$,
the probability for an adsorbed HCO radical to be formed (34.8 \%) is about twice as small as
for CO-(H$_{2}$O)$_{4}$ and CO-(H$_{2}$O)$_{5}$. However, this is also the only species
displaying a significant amount of both Variant 3 and 5, 32.4 \%~and 20.8 \%, respectively.
Furthermore, the likelihood of Variant 3 and 1 are equivalent as 32.4 \%~of the trajectories lead
to an H atom stuck on the water molecules.  This results from the larger surface area
of CO-(H$_{2}$O)$_{10}$ as compared to the other aggregates. It should also be noted that
in those cases, the hydrogen does not diffuse to the CO molecule in the time length of the simulations.

In order to get more dynamical, \textit{i.e.} time-dependent, insights into the H + CO
recombination mechanism, we calculated the probability for HCO formation along the MD simulations. 
Two different functions, P1(t) and P2(t), were computed from all the trajectories, regardless of 
the corresponding Variants. They are defined as follow: -(i)- P1(t) is the cumulative probability that
the HCO radical was formed for the first time at time \textit{t} during the simulation without considering
any further dissociation or recombination; -(ii)- P2(t) is the cumulative probability that HCO was
formed a time \textit{t} and remains stable till the end of the simulation. From these definitions,
P1(t) is mainly influenced by Variants 1, 2 and 7, at a lesser extent it is also influenced by the other
Variants, and P2(t) by trajectories belonging to Variants 1 and 2, only. Figure~\ref{time_func} displays
the P1(t) and P2(t)  functions for CO-(H$_{2}$O)$_{n}$ (n=1, 3, 5 and 10).

\begin{figure}[h!]
\includegraphics[scale=1.2]{new_time.png}
\caption{Cumulative probability of HCO formation as a function of time for the CO-(H$_{2}$O)$_{n}$
(n=1, 3, 5 and 10) clusters.} \label{time_func}
\end{figure}

We can see different dynamical pictures depending on the species.  For CO-(H$_{2}$O),
formation of the HCO radical occurs very rapidly. The curve P1(t)  reaches its maximum
during the first two picoseconds of simulation. Then, it retains a constant  value of 11.6 \%.
This latter value correlates with the sum of the probabilities of Variants 2 and 7  in the
Table~\ref{result_ratio} (10.6 \%). It should be noted that the small difference between
these values results from the definition of both P1(t) and the Variants. Indeed, P1(t) does
not take into account the possible dissociation of HCO and subsequent disconnection of
the reactant from the water cluster. Those mechanisms correspond to Variants 4 to 6. Thus,
in 1\% \ of the trajectories, which belong to Variants 4, 5 or 6, there is the transient
formation of HCO before dissociation and desorption of one or both reactants.
P2(t) shows that a stable radical is formed only in 4.5 \% of cases.
However, this function increases only during the last picosecond of the trajectory. This
demonstrates that the HCO-(H$_{2}$O) aggregate is in a regime of successive formation-dissociation.
The increase of P2(t) at the end of the simulation thus results from the finite time length of the simulations.

For CO-(H$_{2}$O)$_{3}$ and CO-(H$_{2}$O)$_{5}$, other behaviors are observed. Indeed, in the former case,
HCO formation occurs gradually from 0.5 to 6~ps. This strongly contrast with CO-(H$_{2}$O). The maximum
value of P1(t) (66.5 \%) corresponds to the sum of the probabilities of Variants 1, 2 and 7, which is equal to 66.4 \%.
Consequently, a negligible amount of trajectories where HCO is formed lead to its dissociation and the desorption
of the reactants. In addition, in contrast with CO-(H$_{2}$O), a large amount of the successfully formed HCO 
remain stable
until the end of the trajectories. 61.5 \%~of them, corresponding to Variants 1 and 2,  lead to a stable radical. Again,
the increase of the P2(t) curve during the last picosecond of simulation is only an artefact.
In CO-(H$_{2}$O)$_{5}$, a faster HCO formation is observed as revealed by the sharp rise of the P1(t)
curve between 0.5 and 1~ps. After 1~ps, it smoothly increases till 8~ps. The probability of HCO  formation is
63 \%, which correlates with the corresponding Variants in the Table~\ref{result_ratio} (V.1+V.2+V.7 = 62.2 \%).
Similarly, the maximum of the P2(t) curve, 62.1 \%, correspond to the sum of the Variants 1 and 2 (61.4 \%).
Consequently, most of the radicals that are formed remain stable during the considered time length. 
The curves for CO-(H$_{2}$O)$_{3}$ and CO-(H$_{2}$O)$_{5}$ in the Figure~\ref{time_func} display a
time lag between the increase of the P1(t) and P2(t) functions. The faster rise of P1(t) compared to P2(t)
shows that HCO is formed in a number of trajectories but immediately dissociates due to the excess
kinetic energy of the proton. This reveals that a small amount of time is needed to dissipate this energy,
which then allows for the formation of stable HCO and the rise of P2(t). This rise is faster for
CO-(H$_{2}$O)$_{5}$ due to the larger number of water molecules contributing to energy dissipation.
Finally, a plateau is observed for the P1(t) curve of CO-(H$_{2}$O)$_{5}$  at $\sim$1.2~ps. The subsequent
increase of P1(t) can be interpreted as the reaction of CO with hydrogen atoms initially stuck to
water molecules. A visual analysis of the trajectories reveal that only hydrogen atoms that are initially
very close to CO can lead to HCO formation as no diffusion over several water molecules is observed.
This is facilitated by the re-orientation of the CO molecule that allows for the hydrogen capture.
However, this process is rather long which explains the slow rise of the P1(t) and P2(t) curves for
CO-(H$_{2}$O)$_{5}$. The CO-(H$_{2}$O)$_{10}$ curves are very similar to the ones of CO-(H$_{2}$O)$_{5}$
which suggests a similar behavior. The two main differences between those two species are: 
a small time lag for the rise of P1(t) for CO-(H$_{2}$O)$_{10}$, which is attributed to a different initial
orientation of CO with respect to the colliding hydrogen, and a lower intensity of P1(t) and P2(t), attributed
to the larger amount of trajectories leading to the sticking of the hydrogen on the surface.

Finally, we determined the  reaction cross section for the investigated H + CO-(H$_{2}$O)$_{n}$ 
(n=0-5, 10) reaction, defined as:

\begin{eqnarray}\label{cross}
\Omega=\int_{0}^{b_{max}} p(b)2\pi bdb
\end{eqnarray}
where $b$ is the impact parameter and $p(b)$ the probability for a successful HCO formation at
impact parameter $b$. We considered as successful simulations corresponding to Variants 1 and 2,
only. Figure~\ref{size} shows the evolution of the reaction cross section $\Omega$ with increasing
cluster size from 0 to 10 water molecules.

\begin{figure}[h!]
\includegraphics[scale=1.0]{cs_size.png}
\caption{Reaction cross section for H + CO formation as a function of  water cluster size.} \label{size}
\end{figure}

In the case of a single CO molecule, the reaction cross section is null, since a stable HCO radical 
is never formed. When we go from one to two water molecules, although the size of CO-H$_{2}$O 
and CO-(H$_{2}$O)$_{2}$ is similar (the distance between the two most distant atoms is 4.22  and
4.19 \AA, respectively), a significant increase of $\Omega$ is observed. Again, this shows that
the second water molecule assists the capture of the incident hydrogen and makes the
collision-induced reaction more effective. These data are consistent with the previous
analysis of the trajectories.
From two to four  water molecules, the reaction cross section does not significantly change. For five
water  molecules, a second small increase of $\Omega$ can be observed. Consequently, up to 5
water molecules,  the present results suggest that the reaction cross section increases with cluster
size, and hence, the probability of the HCO radical formation also increases.
In the case of CO-(H$_{2}$O)$_{10}$, as previously mentioned for the probabilities of
the collisional pathways, the surface area of (H$_{2}$O)$_{10}$ significantly increases and the probability 
for a successful reaction between H and CO drastically decreases. Indeed, despite considering
frontal collisions only and no orientational sampling, $\Omega$ decreases to $\sim$3~\AA.
This precludes from calculating a cross section to be directly compared with the smaller clusters.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions} \label{Concl}

In the present study, a statistical analysis of the collision-induced reaction between hydrogen 
and carbon monoxide on water clusters (H$_{2}$O)$_{n}$ (n=0-5,10) was achieved 
using molecular dynamics simulation in combination with SCC-DFTB calculations of the energies
and forces. SCC-DFTB allows to describe covalent bonding in a quantum mechanical scheme
and charge relaxation along the dynamics and is thus well suited to study reactivity.
We show that the SCC-DFTB approach allows to satisfactorily describe the interaction between
the  weakly bounded CO and H$_{2}$O molecules and can further be used to reaction study.
In a subsequent step, the classical collision-induced reaction of H  with CO in CO-(H$_{2}$O)$_{n}$
was investigated dynamically and the various pathways for the collision were analysed statistically.
We show that stable HCO radicals can be formed  in species containing two or more water molecules.
Those water molecules play a key role in the collisional mechanism as they allow to dissipate the excess
kinetic energy of the colliding hydrogen.
However, in the case of (H$_{2}$O)$_{2}$ and (H$_{2}$O)$_{3}$, HCO returns to the gas phase in most cases
as dissociation between the radical and the water molecules is  almost always observed. However,
starting from four water molecules, the successfully formed HCO remains associated 
with the surface. This is an important feature for further hydrogenation steps that would lead to methanol.
For three water molecules and above, a large number of collisions lead  to the sticking of the hydrogen
atom on the water molecules. In most cases, no subsequent HCO formation is observed in the simulation 
time length. This suggests that the diffusion of H on the water cluster and its subsequent meeting with CO
(referred to as the Langmuir-Hinshelwood mechanism) either is damped by the cluster character of the
surface (with respect to perfect ice) or occurs at longer time than described by the present simulations.  
In a few cases, when H is initially stuck close to CO, recombination can occur thanks to the
re-orientation of the CO molecule on the cluster allowing for the subsequent formation of HCO. This
process is observed mainly for (H$_{2}$O)$_{5}$ and (H$_{2}$O)$_{10}$. The recombination cross section
is of the order of  $\approx$ 4 to 5~\AA \ for systems containing 3-5 water molecules. When increasing
the number of water molecules up to ten, this value drops drastically due to the sticking of H to the
cluster. Finally, it is worth point out that the size of the considered aggregates as well as the initial
conditions and length time of the trajectories favour the Eley-Rideal process with respect to the HCO
formation following the Langmuir-Hinshelwood mechanism, \textit{i.e.} the diffusion of hydrogen over the
water cluster. This latter would require a completely different set-up of the simulations to be properly
described.

Certainly, the present quantitative results are subject to the accuracy of the SCC-DFTB method and the 
limitation of the dynamical parameters, such as the number of trajectories and their time length. 
Nevertheless, we think that the present results are statistically meaningful and provide clear information
on the influence of water molecules on the H + CO reaction process in a dynamical picture.
The present work 
could be extended in the following directions: -(i)- increase the CO coverage of the water clusters, 
which would avoid the isolation of CO  on a large nanodroplet, possibly with a factor dependence 
of the surface area of the nanodroplet; -(ii)- perform simulations on real ice surface, ordered or
disordered, using periodic boundary conditions.  We hope that the present work brings some
confirmation and new insights for the hydrogenation chain of CO up to methanol and contribute to
the understanding of the chemistry of the ISM.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The "Acknowledgement" section can be given in all manuscript
%% classes.  This should be given within the "acknowledgement"
%% environment, which will make the correct section or running title.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgement}
The authors acknowledge the supercomputing facility of CALMIP for generous allocation
of computer resources (projects P1320 and P0059). The authors declare that there has
been no significant financial support for this work.
\end{acknowledgement}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The same is true for Supporting Information, which should use the
%% suppinfo environment.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% The appropriate \bibliography command should be placed here.
%% Notice that the class file automatically sets \bibliographystyle
%% and also names the section correctly.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{Ref}

\end{document}