These_linjie_JC/thesis/3/structure_stability.tex
2021-06-15 08:35:07 +02:00

723 lines
105 KiB
TeX
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

%Conclusions% this file is c\emph{•}alled up by thesis.tex
% content in this file will be fed into the main document
%%\chapter{Aims of the project} % top level followed by section, subsection
\raggedbottom
\ifpdf
\graphicspath{{3/figures/PNG/}{3/figures/PDF/}{3/figures/}}
\else
\graphicspath{{3/figures/EPS/}{3/figures/}}
\fi
\chapter{Exploration of Structural and Energetic Properties} \label{chap:structure}
This \textbf{third chapter} of my thesis merges two independent studies dealing with the determination of the low-energy isomers of
ammonium/ammonia water clusters, (H$_2$O)$_{n}${NH$_4$}$^+$ and (H$_2$O)$_{n}$NH$_3$, and protonated uracil water clusters, (H$_2$O)$_{n}$UH$^+$.
As highlighted in the general introduction of this thesis and in chapter~\ref{chap:comput_method}, performing global optimization of
molecular clusters is not straightforward. The two studies presented in this chapter thus share a main common methodology which is the
combination of the \textbf{self-consistent-charge density functional based tight-binding} (SCC-DFTB) method for the efficient calculation of the potential
energy surfaces (PES) and the \textbf{parallel-tempering molecular dynamics} (PTMD) approach for their exploration. All low-energy isomers
reported in this chapter are discussed in terms of structure, relative energy and binding energy which are compared to the literature
when available. Calculations at higher level of theory are also performed to refine the results obtained at the SCC-DFTB level or to
validate the results it provides. In particular, in this chapter, an improve set of parameters is proposed to describe sp$^3$ nitrogen
containing compounds at the SCC-DFTB level. Our results are also used to complement collision-induced dissociation experiments
performed by S. Zamith and J.-M. L'Hermite at the \textit{Laboratoire Collisions Agr\'egats R\'eactivit\'e} (LCAR).
%The background of the systems studied in this thesis and theoretical methods related to the projects of this thesis have been introduced in chapter \ref{chap:general_intro} and chapter \ref{chap:comput_method}. Then the two kinds of projects investigated in my Phd will be shown in the following chapter \ref{chap:structure} and chpater \ref{chap:collision}.
%For the topics about ammonium/ammonia water clusters and protonated uracil water clusters, both of the two studies are about the structural and energetic properties. For instance, some lowest-energy isomers were obtained and the relative binding energies of the lowest-energy isomers were calculated. So the work of these two topics are put in the same chapter. These two studies demonstrates the accuracy of SCC-DFTB to describe ammonium/ammonia water clusters and protonated uracil water clusters. The calculation methods for the topics about ammonium/ammonia water clusters and protonated uracil water clusters are almost the same, therefore, the methods of these two topics will be described together. But their introduction, results, and conclusion will be shown independently (in different subsections). For the topic about ammonium/ammonia water clusters, it is a pure theoretical calculation while for the topic about protonated uracil water clusters, it is collaborated with the experiment and some experimental results will be displayed.
\section{Computational Details} \label{sec:structure-methods}
\subsection{SCC-DFTB Potential}
SCC-DFTB electronic structure calculations presented in this chapter were all performed with the deMonNano code.\cite{deMonNano2009} The details
of the method are presented in section~\ref{sec:DFTB} of chapter~\ref{chap:comput_method}.
%DFTB is an approximated DFT scheme whose computational efficiency depends on the use of parameterized integrals.\cite{Elstner2014, Elstner1998, dftb1, dftb2}
The mio-set for the Slater-Koster tables of integrals was used.\cite{Elstner1998} However, it has been shown that these integrals do not properly described
sp$^3$ hybridized nitrogen, in particular, proton affinity.\cite{Gaus2013para} Consequently, in order to avoid spurious deprotonation of the sp$^3$ hybridized
nitrogen in NH$_4^+$ and to correctly reproduce binding energies calculated at the MP2/Def2TZVP level, I propose to modify the original mio-set for
Slater-Koster tables of N-H integrals by applying them a multiplying factor. Several of them were tested and I present here the results obtained for two of them: 1.16 and 1.28. For the sp$^2$ nitrogen of uracil, the original integrals of the mio-set were used. To improve
description of the intermolecular interactions, the original Mulliken charges were replaced by the CM3 charges,\cite{Li1998, Thompson2003, Rapacioli2009corr}
(see equation~\ref{CM3} in section~\ref{sec:DFTB}) and an empirical correction term (see equation~\ref{dispersionE} in section~\ref{sec:DFTB}) was used
to describe dispersion interactions.\cite{Rapacioli2009corr, Elstner2001, Zhechkov2005} A. Simon \textit{et al.} developed a SCC-DFTB potential that leads to geometries, frequencies, and relative energies close to the corresponding experimental and CCSD(T)/aug-cc-pVTZ results.\cite{Simon2012, Odutola1980} The corresponding $D_{OH}$ parameter, \textit{i.e.} 0.129, is retained in the studies presented in thus chapter. $D_{NH}$ is tested in the study of ammonium/ammonia water clusters and two values were retained and thoroughly tested: 0.12 and 0.14. D$_{NO}$ is set to zero.
\subsection{SCC-DFTB Exploration of PES}
To determine the lowest-energy isomers of (H$_2$O)$_{1-10,20}${NH$_4$}$^+$, (H$_2$O)$_{1-10}$NH$_3$
and \newline (H$_2$O)$_{1-7,11,12}$UH$^+$ clusters, their PES were thoroughly explored using PTMD \cite{Sugita1999, Sugita2000, Earl2005}
simulations combined with a SCC-DFTB \cite{Elstner1998} description of the energies and gradients. I describe
below the \textbf{detailed parameters used for all the simulations conducted within this chapter}.
\textbf{Detailed parameters for PTMD simulations of (H$_2$O)$_{1-10,20}${NH$_4$}$^+$ and (H$_2$O)$_{1-10}$NH$_3$ clusters} are as follows.
For (H$_2$O)$_{1-3}${NH$_4$}$^+$ and (H$_2$O)$_{1-3}${NH$_3$} clusters, 16 replicas were used
with a linear distribution of temperatures with a 15~K step ranging from 10 to 250 K. 40 replicas with a 6~K step ranging from 10 to 250~K were
considered for (H$_2$O)$_{4-10,20}${NH$_4$}$^+$ and (H$_2$O)$_{4-10}${NH$_3$} species. All trajectories were 5~ns long and a time step
of 0.5 fs was used to integrate the equations of motion. A Nos{\'e}-Hoover chain of 5 thermostats was employed for all the simulations to achieve
simulations in the canonical ensemble.\cite{Nose1984M, Hoover1985} Thermostat frequencies were fixed at 400 cm$^{-1}$.
To identify low-energy isomers of (H$_2$O)$_{1-3}${NH$_4$}$^+$ and (H$_2$O)$_{1-3}${NH$_3$} clusters, 303 geometries were periodically
selected from each replicas and further optimized at the SCC-DFTB level, which produced 4848 optimized geometries per cluster. For
(H$_2$O)$_{4-10,20}${NH$_4$}$^+$ and (H$_2$O)$_{4-10}${NH$_3$} clusters, 500 geometries were periodically selected from each
replicas leading to 20000 optimized geometries per cluster. For (H$_2$O)$_{20}${NH$_4$}$^+$, the initial structure used for the global
optimization process was the lowest-energy structure reported by J. Douady \textit{et al.}.\cite{Douady2009} The five lowest-energy
isomers among the 4848 or 20000 optimized geometries were further optimized using the MP2/Def2TZVP method. See below for the
details on MP2/Def2TZVP calculations.
\textbf{Detailed parameters for PTMD simulations of (H$_2$O)$_{1-7,11,12}$UH$^+$ clusters} are as follows.
40 replicas with temperatures rnaging linearly from 50 to 350 K were used. Each trajectory was 4 ns long, and the integration time step was 0.5 fs.
A reasonable time interval for the PT exchanges was 2.5 ps. A Nos{\'e}-Hoover chain of five thermostats with frequencies of 800 cm$^{-1}$ was
applied to achieve an exploration in the canonical ensemble.\cite{Nose1984M, Hoover1985} To avoid any spurious influence of the initial
geometry on the PES exploration, three distinct PTMD simulations were carried out with distinct initial proton location: on the uracil in two cases
and on a water molecule in the other one. In the former cases, I used two isomers u178 and u138 of UH$^+$ shown in Figure~\ref{uracil_s} as
the initial geometries.\cite{Wolken2000, Pedersen2014} 600 geometries per temperature were linearly selected along each PTMD simulation
for subsequent geometry optimization leading to 72000 structures optimized at SCC-DFTB level. These structures were sorted in ascending
energy order and checked for redundancy. 9, 23, 46, 31, 38, 45, 63, 20, and 29 structures were then selected for (H$_2$O)UH$^+$,
(H$_2$O)$_{2}$UH$^+$, (H$_2$O)$_{3}$UH$^+$, (H$_2$O)$_{4}$UH$^+$, (H$_2$O)$_{5}$UH$^+$, (H$_2$O)$_{6}$UH$^+$,
(H$_2$O)$_{7}$UH$^+$, (H$_2$O)$_{11}$UH$^+$ and (H$_2$O)$_{12}$UH$^+$ respectively, to perform geometry optimizations at the MP2/Def2TZVP
level. See below for the details on MP2/Def2TZVP calculations.
%\figuremacro{uracil}{Structures of the two protonated uracil isomers, u178 (keto-enol form) and u138 (di-keto form), used as initial conditions in the PTMD simulations.}
\begin{figure}[h!]
\includegraphics[width=0.5\linewidth]{uracil.pdf}
\centering
\caption{Structures of the two protonated uracil isomers, u178 (keto-enol form) and u138 (di-keto form), used as initial conditions in the PTMD simulations.}
\label{uracil_s}
\end{figure}
\subsection{MP2 Geometry Optimizations, Relative and Binding Energies}
Some low-energy isomers obtained at the SCC-DFTB level were further optimized at the MP2 level of theory in combinationwith an all electron Def2TZVP basis-set.\cite{Weigend2005, Weigend2006} All calculations used a a tight criteria for geometry
convergence and an ultrafine grid for the numerical integration. All MP2 calculations were performed with the Gaussian 09
package.\cite{GaussianCode}
\textbf{Detailed parameters for (H$_2$O)$_{1-10,20}${NH$_4$}$^+$ and (H$_2$O)$_{1-10}$NH$_3$ clusters.}
Following SCC-DFTB optimizations, the five lowest-energy isomers of (H$_2$O)$_{1-10}${NH$_4$}$^+$ and
(H$_2$O)$_{1-10}${NH$_3$} clusters were further optimized at the MP2/Def2TZVP level of theory.
In section~\ref{sec:ammoniumwater}, relatives energies with respect to the lowest-energy isomer of each
cluster are reported. Impact of zero-point vibrational energy (ZPVE) corrections on relative
energies were evaluated at MP2/Def2TZVP level. To evaluate the strength of water-ammonium and water-ammonia
interactions and to assess the accuracy of the SCC-DFTB method, binding energies are also reported.
Two distinct approaches were used to calculate binding energies. The first one considers only the binding
energy between the water cluster as a whole and the impurity, {NH$_4$}$^+$ or NH$_3$, while the second one
considers the binding energy between all the molecules of the cluster. In both cases, the geometry of the molecules
is the one found in the optimized cluster. Using these two methods, relative binding energies (E$_{bind.}$(SCC-DFTB)-E$_{bind.}$(MP2/Def2TZVP))
$\Delta E_{bind.}^{whole}$ and $\Delta E_{bind.}^{sep.}$ were obtained. For all binding energies of
(H$_2$O)$_{1-10}${NH$_4$}$^+$ and (H$_2$O)$_{1-10}${NH$_3$} clusters calculated at MP2/Def2TZVP level,
basis set superposition errors (BSSE) correction was considered by using the counterpoise method of Boys and
Bernardi.\cite{Boys2002}
\textbf{Detailed parameters for (H$_2$O)$_{1-7, 11, 12}$UH$^+$ clusters.}
Following SCC-DFTB optimizations, the six lowest-energy isomers of
(H$_2$O)$_{1-7, 11, 12}$UH$^+$ clusters were further optimized at the MP2/Def2TZVP level of theory.
The binding and relative energies calculated at MP2/Def2TZVP level without BSSE correction
of clusters (H$_2$O)$_{2-7, 11, 12}$UH$^+$ are discussed in section~\ref{structureUH}.
\subsection{Structure Classification}
For clusters (H$_2$O)$_{1-10}${NH$_4$}$^+$ and (H$_2$O)$_{1-10}${NH$_3$},
\textquotedblleft n-x\textquotedblright and \textquotedblleft n$^\prime$-x\textquotedblright labels are used to distinguish between the
(H$_2$O)$_{n}${NH$_4$}$^+$ and (H$_2$O)$_{n}${NH$_3$} reported isomers, respectively, obtained at the SCC-DFTB level. For comparison,
\textquotedblleft n-x\textquotedblright and \textquotedblleft n$^\prime$-x\textquotedblright isomers are also optimized at the MP2/Def2TZVP
level. In that case, the resulting structures are referred to as \textquotedblleft n-x$^*$\textquotedblright and
\textquotedblleft n$^\prime$-x$^*$\textquotedblright~ to distinguish them more easily although they display the same general topology as
\textquotedblleft n-x\textquotedblright and \textquotedblleft n$^\prime$-x\textquotedblright isomers.
In these notations, n and n$^\prime$ denote the number of water molecules in the ammonium and ammonia water clusters, respectively.
x is an alphabetic character going from a to e that differentiates between the five low-energy isomers reported for each cluster in ascending
energy order, \textit{i.e.} a designates the lowest-energy isomer.
%\section{\label{sec:ammoniumwater}The structure and energetics properties study of ammonium or ammonia including water clusters}
\section{Structural and Energetic Properties of Ammonium/Ammonia including Water Clusters} \label{sec:ammoniumwater}
\subsection{General introduction}
%Ionic clusters (typically made up of a core ion surrounded by one or more solvating molecules) are known to be involved in the chemistry of the upper and mid atmosphere [6].《Ion rearrangement at the beginning of cluster formation: isomerization pathways and dissociation kinetics for the ionized dimethylamine dimer》
Water clusters play an important role in various areas such as atmospheric and astrochemical science, chemistry and biology.\cite{Keesee1989, Gilligan2000,
Sennikov2005, Cabellos2016, Orabi2013, Bommer2016, Rodgers2003, Van2004, Gibb2004, Tielens2005, Parise2005, Boogert2015, Dulieu2010, Michoulier2018}
They are involved into the critical stages of nucleation and growth of water-containing droplets in the atmosphere thus contributing to the physical and
chemical properties of this medium.\cite{Kulmala2004}
%For instance, water clusters are able to absorb a significant amount of radiant energy thus decreasing the greenhouse effect.\cite{Galashev2010}
In many cases, the presence of chemical impurities interacting with water aggregates strongly affect their properties.
For instance, ammonia is an important compound commonly found in the atmosphere and which displays a key role in aerosol chemistry.\cite{Ziereis1986} Its high
basicity makes it a potential proton sink that can form a ionic center for nucleation.\cite{Perkins1984, Arnold1997} E. Dunne and co-workers also reported that most nucleation occurring in the atmosphere involves ammonia or biogenic organic compounds, in addition to sulfuric acid.\cite{Dunne2016} J. Kirkby \textit{et al.}
also found that even a small amount of atmospherically relevant ammonia can increase the nucleation rate of sulphuric acid particles by several orders of magnitude.\cite{Kirkby2011}
The significance of ammonium and ammonia water clusters have thus motivated a large amount of experimental and theoretical studies during the past decades.\cite{Perkins1984, Herbine1985, Stockman1992, Hulthe1997, Wang1998, Chang1998, Jiang1999, Hvelplund2010, Douady2009, Douady2008, Morrell2010, Bacelo2002, Galashev2013}
As a few examples, in 1984, (H$_2$O)$_{2}${NH$_4$}$^+$ was identified using mass spectrometry by M. D. Perkin \textit{et al.}\cite{Perkins1984} In 1997, Stenhagen
and co-workers studied the {(H$_2$O)$_{20}$H$_3$O}$^+$ and (H$_2$O)$_{20}${NH$_4$}$^+$ clusters and found that both species display a similar
structure.\cite{Hulthe1997} P. Hvelplund \textit{et al.} later reported a combined experimental and theoretical study devoted to protonated mixed ammonia/water
which highlighted the idea that small protonated mixed clusters of water and ammonia contain a central {NH$_4$}$^+$ core.\cite{Hvelplund2010}
%The (H$_2$O)NH$_3$ complex has been experimentally investigated via radio frequency and microwave spectra by Herbine \textit{et al.}, and via microwave and tunable far-infrared spectroscopy by Stockman and co-workers.\cite{Herbine1985, Stockman1992}
Theoretical calculations devoted to ammonium and ammonia water clusters have also been extensively conducted.\cite{Lee1996, Chang1998,
Skurski1998, Jiang1999, Donaldson1999, Sadlej1999, Hvelplund2010, Bacelo2002, Galashev2013} Among them, J. Novoa \textit{et al.} studied the (H$_2$O)$_4$NH$_3$
aggregate and found the existence of a minimum in its potential energy surface corresponding to a (H$_2$O)$_{3}$···{NH$_4$}$^+$···OH$^-$ structure, resulting from one proton transfer from a water molecule to the ammonia molecule.\cite{Lee1996} D. Bacelo later reported a number of low-energy minima for (H$_2$O)$_{3-4}$NH$_3$ clusters obtained from \textit{ab initio} calculation and a Monte Carlo exploration of the PES.\cite{Bacelo2002} More recently, J. Douady \textit{et al.} performed a global optimization of (H$_2$O)$_{n}${NH$_4$}$^+$ ($n$ = 1-24) clusters again using a Monte Carlo procedure in combination with a Kozack and Jordan empirical force field.\cite{Douady2008, Kozack1992polar}
In this study, the finite temperature properties as well as vibrational signature of several clusters thus highlighting the key contribution of simulations in understanding such species. Morrell and Shields also studied the
(H$_2$O)$_{n}${NH$_4$}$^+$ ($n$ = 1-10) aggregates \textit{via} a mixed molecular dynamics and quantum mechanics methodology to calculate energies and free energies
of formations which were in good agreement with previous experimental and theoretical results.\cite{Morrell2010}
More recently, S. Pei \textit{et al.} determined that (H$_2$O)$_{n}$NH$_4^+$ clusters start to adopt a closed-cage geometry at $n$=8.\cite{Pei2015}
Finally, W. Walters and collaborators determined the geometry of (H$_2$O)$_{16}$NH$_3$ and (H$_2$O)$_{16}$NH$_4^+$ at the HF/6-31G(d) level,
and observed strong hydrogen bonding between water and the lone pair of NH$_3$ and bewteen NH$_4^+$ and the four adjacent water molecules.\cite{Walters2018}
As for the study of other molecular clusters, the range of applicability of theoretical simulations to describe ammonium and ammonia water clusters is dictated
by the balance between accuracy, transferability and computational efficiency. While \textit{ab-initio} methods can accurately model small aggregates, their
application to large species is more difficult, in particular when an exhaustive exploration of the PES is required. In contrast, force-field
potentials are computationally extremely efficient and can be coupled to global optimization methods but their transferability is limited.
The SCC-DFTB approach can be seen as an intermediate approach which combines the strengths of both \textit{ab-initio} and force-field methods.
Indeed, it can be as accurate as DFT while computationally more efficient and is more transferable than force fields (see chapter \ref{chap:comput_method})
In recent years, SCC-DFTB has been successfully applied to the study of various molecular clusters: pure, protonated, and de-protonated water
clusters,\cite{Choi2010, Choi2013, Korchagina2017, Simon2019} water clusters on PAHs,\cite{Simon2012, Simon2013water}
sulfate-containing water clusters,\cite{Korchagina2016} water clusters in an argon matrix,\cite{Simon2017formation} whether it is for global optimization or
for the study of finite-temperature properties. However, in its original formulation, SCC-DFTB does not provide good results for the description of ammonia
and ammonium as nitrogen hybridization seems to be a problem for minimal basis-set methods like SCC-DFTB.\cite{Winget2003} M. Elstner and coworkers found
consistent errors (about 14.0 kcal.mol$^{-1}$) for deprotonation energies of sp$^3$ hybridized nitrogen containing systems, whereas sp$^1$ and sp$^2$ systems
display much smaller errors.\cite{Gaus2013para}
In this section, I first propose an improvement of the SCC-DFTB scheme to describe ammonium and ammonia water clusters by modifying
both Hamiltonian and overlap N-H integrals and introducing optimized atomic charges.\cite{Thompson2003, Rapacioli2009} By combining this
improved SCC-DFTB scheme with PTMD simulations, global optimization of the (H$_2$O)$_{1-10}${NH$_4$}$^+$ and (H$_2$O)$_{1-10}${NH$_3$}
clusters is then performed which allows to report a number of low-energy isomers for these species. Among them, a selected number of structures
are further optimized at the MP2/Def2TZVP level of theory to confirm they are low-energy structures of the PES and to rationalize the difference in
relative energy between both methods. A detailed description of the reported low-energy isomers is then provided as well as comparisons with the
literature. The heat capacity curve of (H$_2$O)$_{20}${NH$_4$}$^+$ is also obtained at the SCC-DFTB level and compared to previously published
simulations. Some conclusions are finally presented. A very small part of this work has been published in 2019 in a review in\textit{Molecular Simulation}.\cite{Simon2019}
A full paper devoted to this work is in preparation.
\subsection{Results and Discussion}
\subsubsection{Dissociation Curves and SCC-DFTB Potential}
In order to define the best SCC-DFTB parameter to model ammonia and ammonium water clusters, I have tested various sets of corrections.
Each correction involves two modifications of the potential, the first one is the CM3 charge parameter D$_{NH}$ and the second one is the
multiplying factor, noted $xNH$, applied to the NH integrals in the Slater-Koster tables. So a given set is noted D$_{NH}$/$xNH$. Two sets of
corrections have provided satisfactory results, 0.12/1.16 and 0.14/1.28. Figure~\ref{fig:E_nh4} and ~\ref{fig:E_nh3} present dissociation
curves obtained at the MP2/Def2TZVP, MP2/Def2TZVP with BSSE correction, original SCC-DFTB, SCC-DFTB 0.14/1.28 and SCC-DFTB
0.12/1.16 levels of theory. These curves are obtained using the same set of geometries regardless of the method applied to calculate
the binding energies. They are obtained from the MP2/Def2TZVP optimized structures in which the distance between the water and the
ammonium/ammonia was shifted along the N--O vector, all other geometrical parameters being kept fixed.
\begin{figure}[h!]
\includegraphics[width=0.6\linewidth]{E-distance-nh4-w.png}
\centering
\caption{Binding energies of (H$_2$O){NH$_4$}$^+$ as a function of the N---O distance at MP2/Def2TZVP (plain black), MP2/Def2TZVP with BSSE correction (dotted black),
original SCC-DFTB (plain red), SCC-DFTB (0.14/1.28) (dotted red) and SCC-DFTB (0.12/1.16) (dashed red) levels of theory.}
\label{fig:E_nh4}
\end{figure}
\begin{figure}[h!]
\includegraphics[width=0.6\linewidth]{E-distance-nh3-w.png}
\centering
\caption{Binding energies of (H$_2$O){NH$_3$} as a function of the N---O distance at MP2/Def2TZVP (plain black), MP2/Def2TZVP with BSSE correction (dotted black),
original SCC-DFTB (plain red), SCC-DFTB (0.14/1.28) (dotted red) and SCC-DFTB (0.12/1.16) (dashed red) levels of theory.}
\label{fig:E_nh3}
\end{figure}
From Figure~\ref{fig:E_nh4}, the five curves display the same trends with a minimum located at almost the same N---O distance. At the curve minimum,
binding energies vary between -25.57 and -21,07~kcal.mol$^{-1}$ at the original SCC-DFTB and SCC-DFTB 0.14/1.28 levels, respectively. The binding
energy obtained at the SCC-DFTB 0.12/1.16 level is the closest to that obtained at MP2/Def2TZVP level with BSSE correction with a binding energy difference of
only 0.47~kcal.mol$^{-1}$. The SCC-DFTB 0.14/1.28 curve is also very close with a difference in binding energy only 0.16~kcal.mol$^{-1}$ higher. It is
worth mentioning that both sets of corrections lead to improved results as compared to the original SCC-DFTB parameters which leads to a too low
binding energy as compared to MP2/Def2TZVP level with BSSE correction. Also the position of the minimum is more shifted at the original SCC-DFTB
level (2.64~\AA) than with corrections (2.73~\AA). So from structural and energetic point of views, both sets of corrections are satisfactory.
From Figure~\ref{fig:E_nh3}, the five curves display significant differences. This effect is accentuated by smaller binding energy values: they
vary from -3.82 to -7,39~kcal.mol$^{-1}$ at the original SCC-DFTB and MP2/Def2TZVP levels, respectively, at the minimum of the curves. The binding
energy obtained at the SCC-DFTB 0.12/1.16 level is the closest to that obtained at MP2/Def2TZVP level with BSSE correction with a binding energy
difference of only 0.01~kcal.mol$^{-1}$. The SCC-DFTB 0.14/1.28 curve is also rather close with a difference in binding energy only 1.3~kcal.mol$^{-1}$
higher. Here also, both sets of corrections lead to improved results as compared to the original SCC-DFTB parameters. The position of the minimum
is also well reproduced by the corrected potentials. In contrast to (H$_2$O){NH$_4$}$^+$, the shape of the curves for (H$_2$O){NH$_3$} obtained
at the SCC-DFTB level differs significantly from those obtained at MP2 level. Vibrational frequencies calculated at the SCC-DFTB level for this systems
are therefore expected to be inacurate. It is worth mentioning that the large difference in binding energy between
(H$_2$O){NH$_4$}$^+$ and (H$_2$O)NH$_3$ was expected owing to a stronger electrostatic contribution of {NH$_4$}$^+$ to the
binding energy.
Another very important point when comparing the original SCC-DFTB potential and the corrected potentials, is the structure obtained for the
(H$_2$O){NH$_4$}$^+$ dimer. Figure~\ref{dimers} compares the structure obtained from geometry optimization at the SCC-DFTB 0.14/1.28
and original SCC-DFTB levels. The N-H covalent bond involved in the hydrogen bond is longer with the original potential while the N---O
distance is smaller by 0.14~\AA. This is reminiscent of the too low proton affinity of {NH$_4$}$^+$ predicted by the original SCC-DFTB potential.
This discrepancy has been previously highlighted in other studies,\cite{Winget2003,Gaus2013para} and makes this potential unusable in any
realistic molecular dynamics simulation as it leads to a spurious deprotonation. Both sets of corrections are free of this error.
%\red{An additional
%proof of this assertion is presented in Figure~\ref{XX} that displays the energy barrier for the proton transfer at constant N---O distance obtained
%from different methods.}
\begin{figure}[h!]
\includegraphics[width=0.3\linewidth]{dimers.png}
\centering
\caption{Structure of (H$_2$O){NH$_4$}$^+$ obtained from geometry optimization at the SCC-DFTB 0.14/1.28 (right) and original SCC-DFTB (left) levels.}
\label{dimers}
\end{figure}
Figures~\ref{fig:E_nh4} and ~\ref{fig:E_nh3} show that SCC-DFTB 0.12/1.16 better describe both (H$_2$O){NH$_3$} and
(H$_2$O){NH$_4$}$^+$ dissociation curves. Furthermore, as (H$_2$O){NH$_3$} is characterized by a much lower
binding energy than (H$_2$O){NH$_4$}$^+$, an error of the order of $\sim$1.0~kcal.mol$^{-1}$ is more likely to play
a significant role for ammonia than ammonium containing species. All the following discussion therefore involve the
SCC-DFTB 0.12/1.16 potential.
\subsubsection{Small Species: (H$_2$O)$_{1-3}${NH$_4$}$^+$ and (H$_2$O)$_{1-3}${NH$_3$}}
As a first test case for the application of the SCC-DFTB 0.14/1.28 potential is the study of small ammonium and ammonia water clusters:
(H$_2$O)$_{1-3}${NH$_4$}$^+$ and (H$_2$O)$_{1-3}${NH$_3$}. Due to the limited number of low-energy isomers for these species, we
only consider the lowest-energy isomer of (H$_2$O)$_{1-2}${NH$_4$}$^+$ and (H$_2$O)$_{1-3}${NH$_3$} and the two lowest-energy
isomers for (H$_2$O)$_{3}${NH$_4$}$^+$.
As displayed in Figure~\ref{fig:nh3-nh4-1w} and\ref{fig:nh3-nh4-2-3w}, the reported low-energy isomers 1-a, 1$^\prime$-a, 2-a, 2$^\prime$-a, 3-a, 3-b, and
3$^\prime$ display a structure very similar to those obtained at the MP2/Def2TZVP level (1-a$^*$, 1$^\prime$-a$^*$, 2-a$^*$, 2$^\prime$-a$^*$, 3-a$^*$,
3-b$^*$ and 3$^\prime$-a$^*$). Indeed, although differences in bond lengths are observed, they are rather small.
In terms of energetics,
From en energetic point of view, it is interesting to first look at the relative energy between the two reported isomers of
(H$_2$O)$_{3}${NH$_4$}$^+$. Isomer 3-b is 2.12~kcal·mol\textsuperscript{-1} higher than 3-a at the SCC-DFTB level.
At the MP2/Def2TZVP level, 3-b is 0.30~kcal·mol\textsuperscript{-1} lower than 3-a when ZPVE is not considered while
it is 1.21 kcal·mol\textsuperscript{-1} higher when it is considered. In comparison, in the experimental results by
H. Chang and co-workers, 3-a is more stable than 3-b.\cite{Wang1998, Jiang1999} The authors also complemented their
measurements by theoretical calculations that show that at the B3LYP/6-31+G(d) level, 3-a is higher than 3-b but. In
contrast, at the MP2/6-31+G(d) level corrected with ZPVE, the energy of 3-a is lower than that of 3-b while it is inverted
if ZPVE is taken into account.\cite{Wang1998, Jiang1999} Additionally, F. Spiegelman and co-workers, conducted a
global Monte Carlo optimizations with an intermolecular polarizable potential that lead to 3-a as lowest-energy isomer.\cite{Douady2008}
All these results show that for the specific question of lowest-energy isomer of (H$_2$O)$_{3}${NH$_4$}$^+$, SCC-DFTB
has an accuracy close to other \textit{ab initio} methods which confirms its applicability.
\begin{figure}[h!]
\includegraphics[width=0.2\linewidth]{nh3-nh4-1w.png}
\centering
\caption{Structure of 1-a and 1$^\prime$-a isomers obtained at the SCC-DFTB level and corresponding
structures obtained at MP2/Def2TZVP level (1-a$^*$ and 1$^\prime$-a$^*$ isomers). Selected bond
lengths are in \AA.}
\label{fig:nh3-nh4-1w}
\end{figure}
\begin{figure}[h!]
\includegraphics[width=0.5\linewidth]{nh3-nh4-2w.png}
\centering
\caption{Structure of 2-a and 2$^\prime$-a isomers obtained at the SCC-DFTB level and corresponding
structures obtained at MP2/Def2TZVP level (2-a$^*$, 2$^\prime$-a$^*$ isomers). Selected bond
lengths are in \AA.}
\label{fig:nh3-nh4-2-3w}
\end{figure}
\begin{figure}[h!]
\includegraphics[width=0.5\linewidth]{nh3-nh4-3w.png}
\centering
\caption{Structure of 3-a, 3-b and 3$^\prime$-a isomers obtained at the SCC-DFTB level and corresponding
structures obtained at MP2/Def2TZVP level (3-a$^*$, 3-b$^*$ and 3$^\prime$-a$^*$ isomers). Selected bond
lengths are in \AA.}
\label{fig:nh3-nh4-3w}
\end{figure}
\begin{table*}
\begin{center}
\caption{Relative binding energies $\Delta E_{bind.}^{whole}$ and $\Delta E_{bind.}^{sep.}$
of the low-energy isomers of (H$_2$O)$_{1-3}${NH$_4$}$^+$ and (H$_2$O)$_{1-3}${NH$_3$} clusters.
Values are given in kcal.mol$^{-1}$.} \label{reBindE-small}
\begin{tabular}{c|c|c|c|c|c}
\hline
(H$_2$O)$_{n}$NH$_4^+$ & $\Delta E_{bind.}^{whole}$ & $\Delta E_{bind.}^{sep.}$ & (H$_2$O)$_{n}${NH$_3$} & $\Delta E_{bind.}^{whole}$ & $\Delta E_{bind.}^{sep.}$\\
\hline
\rowcolor{lightgray} 1-a & 1.21 & 1.21 & 1$^\prime$-a & -1.17 & -1.17 \\
2-a & 0.82 & 0.91 & 2$^\prime$-a & 0.57 & 0.28 \\
\rowcolor{lightgray} 3-a & -0.25 & 0.11 & 3$^\prime$-a & 0.91 & 0.01 \\
\rowcolor{lightgray} 3-b & 1.21 & -0.15 & - & - & - \\
\hline
\end{tabular}
\end{center}
\end{table*}
As listed in Table~\ref{reBindE-small}, the relative binding energies $\Delta E_{bind.}^{whole}$ or $\Delta E_{bind.}^{sep.}$ of
(H$_2$O){NH$_4$}$^+$ and (H$_2$O){NH$_3$} are 1.21 and -1.17 kcal·mol\textsuperscript{-1}, respectively, which again
highlights that SCC-DFTB is in agreement with MP2/Def2TZVP. For (H$_2$O){NH$_3$}, the negative value
show that MP2/Def2TZVP binding energy is smaller than the SCC-DFTB value. This is inverse to what is shown in Figure~\ref{fig:E_nh3}
and results from structural reorganization after optimization. All other values of Table~\ref{reBindE-small} are equal of smaller than these
values, whether considering $\Delta E_{bind.}^{whole}$ or $\Delta E_{bind.}^{sep.}$, which again demonstrates that the presently proposed
SCC-DFTB potential provides results in line with reference MP2/Def2TZVP calculations.
(H$_2$O)$_{4-10}${NH$_4$}$^+$ clusters have been studied by molecular dynamic and Monte Carlo simulations in combination with
DFT and MP2 approaches although these latter are computationally expensive.\cite{Wang1998, Jiang1999, Douady2008, Lee2004, Douady2009, Morrell2010}
In contrast, to the best of our knowledge, no theoretical calculation about (H$_2$O)$_{5-10}${NH$_3$} clusters have been conducted.
The low computational cost of SCC-DFTB and its seemingly good performances on small clusters provide an appealing opportunity
to thoroughly explore the PES of both large ammonium and ammonia containing water clusters. In the following section, the five lowest-energy
isomers of clusters (H$_2$O)$_{4-10}${NH$_4$}$^+$ and (H$_2$O)$_{4-10}$NH$_3$ are presented and discussed in details.
\subsubsection{Properties of (H$_2$O)$_{4-10}${NH$_4$}$^+$ Clusters}
The five lowest-energy isomers of (H$_2$O)$_{4}${NH$_4$}$^+$ are depicted in Figure~\ref{fig:nh4-4-6w}. 4-a is the lowest-energy isomer obtained
from the global SCC-DFTB optimization and also the lowest-energy configuration after optimization at MP2/Def2TZVP level with ZPVE
corrections. This result is consistent with previous computational studies\cite{Wang1998, Jiang1999, Douady2008, Lee2004, Pickard2005} and
the experimental studies by H. Chang and co-workers.\cite{Chang1998, Wang1998} Isomer 4-a displays four hydrogen bonds around the ionic
center which lead to no dangling N-H bonds. Other isomers of comparable stability are displayed in Figure~\ref{fig:nh4-4-6w}
The energy ordering of 4-a to 4-e at SCC-DFTB level is consistent with that at MP2/Def2TZVP level with ZPVE correction, although they
are slightly higher by$\sim$2.0 kcal.mol$^{-1}$. Isomer 4-c was not reported in H. Changs study,\cite{Jiang1999} and the corresponding
energy ordering of the five lowest-energy isomers was the same as ours which certainly results from the use of a different basis set.
\begin{figure}[h!]
\includegraphics[width=0.9\linewidth]{nh4-4-6w.png}
\centering
\caption{Five lowest-energy isomers of (H$_2$O)$_{4-6}${NH$_4$}$^+$ and corresponding relative energies at MP2/Def2TZVP level with (bold) and without ZPVE (roman)
correction and SCC-DFTB level (italic).
Relative energies are given in kcal·mol\textsuperscript{-1}.}
\label{fig:nh4-4-6w}
\end{figure}
\begin{table*}
\begin{center}
\caption{Relative binding energies $\Delta E_{bind.}^{whole}$ and $\Delta E_{bind.}^{sep.}$ of the five
lowest-energy isomers of (H$_2$O)$_{4-10}${NH$_4$}$^+$ and (H$_2$O)$_{4-10}${NH$_3$}.
Binding energies are given in kcal·mol\textsuperscript{-1}.} \label{reBindE}
\begin{tabular}{c|c|c|c|c|c}
\hline
(H$_2$O)$_{n}$NH$_4^+$ & $\Delta E_{bind.}^{whole}$ & $\Delta E_{bind.}^{sep.}$ & (H$_2$O)$_{n}${NH$_3$} & $\Delta E_{bind.}^{whole}$ & $\Delta E_{bind.}^{sep.}$\\
\hline
\rowcolor{lightgray} 4-a & -1.67 & -0.87 & 4$^\prime$-a & -1.11 & -1.76 \\
\rowcolor{lightgray} 4-b & 0.00 & 0.61 & 4$^\prime$-b & -0.29 & -1.62 \\
\rowcolor{lightgray} 4-c & 0.77 & 0.44 & 4$^\prime$-c & -0.29 & -1.38 \\
\rowcolor{lightgray} 4-d & 0.77 & 0.42 & 4$^\prime$-d & 1.08 & -0.49 \\
\rowcolor{lightgray} 4-e & -4.04 & 0.69 & 4$^\prime$-e & 1.02 & -1.07 \\
5-a & -1.62 & 0.56 & 5$^\prime$-a & 0.82 & -1.78 \\
5-b & 0.72 & 0.48 & 5$^\prime$-b & -0.23 & -2.26 \\
5-c & 0.69 & 0.55 & 5$^\prime$-c & -0.34 & -2.50 \\
5-d & -1.08 & -0.78 & 5$^\prime$-d & -0.59 & -1.84 \\
5-e & -2.08 & 0.88 & 5$^\prime$-e & -0.38 & -2.60 \\
\rowcolor{lightgray} 6-a & -1.71 & -0.38 & 6$^\prime$-a & -0.27 & -3.05 \\
\rowcolor{lightgray} 6-b & -1.14 & -0.76 & 6$^\prime$-b & -0.31 & -3.55 \\
\rowcolor{lightgray} 6-c & -2.06 & 0.27 & 6$^\prime$-c & -1.11 & -4.67 \\
\rowcolor{lightgray} 6-d & -2.90 & -1.06 & 6$^\prime$-d & -0.05 & -4.44 \\
\rowcolor{lightgray} 6-e & -1.18 & -0.60 & 6$^\prime$-e & 0.55 & -1.96 \\
7-a & -2.95 & -0.39 & 7$^\prime$-a & 1.09 & -2.02 \\
7-b & -2.92 & -0.38 & 7$^\prime$-b & -0.02 & -4.07 \\
7-c & -2.17 & 0.09 & 7$^\prime$-c & -0.40 & -4.15 \\
7-d & -1.28 & -1.35 & 7$^\prime$-d & -0.14 & -3.10 \\
7-e & -3.22 & -2.27 & 7$^\prime$-e & -1.11 & -4.32 \\
\rowcolor{lightgray} 8-a & -2.20 & -1.63 & 8$^\prime$-a & -1.12 & -4.41 \\
\rowcolor{lightgray} 8-b & -1.61 & -2.01 & 8$^\prime$-b & -0.10 & -3.04 \\
\rowcolor{lightgray} 8-c & -3.71 & -1.17 & 8$^\prime$-c & -0.41 & -4.46 \\
\rowcolor{lightgray} 8-d & -2.43 & -0.36 & 8$^\prime$-d & 0.20 & -3.68 \\
\rowcolor{lightgray} 8-e & -0.55 & 0.35 & 8$^\prime$-e & -1.28 & -4.75 \\
9-a & -2.02 & -1.39 & 9$^\prime$-a & -0.15 & -4.47 \\
9-b & 0.51 & -0.84 & 9$^\prime$-b & -1.01 & -4.45 \\
9-c & -3.31 & -0.85 & 9$^\prime$-c & -1.04 & -4.42 \\
9-d & -1.58 & -1.78 & 9$^\prime$-d & -1.09 & -5.14 \\
9-e & -2.39 & -0.91 & 9$^\prime$-e & 0.41 & -2.57 \\
\rowcolor{lightgray} 10-a & -2.64 & -1.94 & 10$^\prime$-a & -0.03 & -4.80 \\
\rowcolor{lightgray} 10-b & -5.79 & -4.35 & 10$^\prime$-b & 0.13 & -5.61 \\
\rowcolor{lightgray} 10-c & -1.26 & -2.36 & 10$^\prime$-c & -0.62 & -6.50 \\
\rowcolor{lightgray} 10-d & -1.98 & -1.42 & 10$^\prime$-d & -1.10 & -6.30 \\
\rowcolor{lightgray} 10-e & -7.17 & -1.54 & 10$^\prime$-e & 0.23 & -8.36 \\
\hline
\end{tabular}
\end{center}
\end{table*}
The relative binding energy of SCC-DFTB method to MP2/Def2TZVP method with BSSE correction for isomers 4-a to 4-e are listed in Table \ref{reBindE}. When the four water molecules are considered as a whole part to calculate the binding energy, the relative binding energy of isomers 4-a to 4-e are -1.67, 0.00, 0.77, 0.77 and -4.04 kcal·mol\textsuperscript{-1}. As shown in Table \ref{reBindE}, for isomers 4-a to 4-e, when the four water molecules are separately considered using the geometry in the cluster to calculate the binding energy, the biggest absolute value of the relative binding energy is 0.87 kcal·mol\textsuperscript{-1}. This shows the results of SCC-DFTB are in good agreement with those of MP2/Def2TZVP with BSSE correction for (H$_2$O)$_{4}${NH$_4$}$^+$. From the relative binding energy of (H$_2$O)$_{4}${NH$_4$}$^+$, it indicates that all the water molecules considered as a whole part or separately has an effect on the relative binding energy for the cluster (H$_2$O)$_{4}${NH$_4$}$^+$ and the overall $\Delta E_{bind.}^{whole}$ are bigger than $\Delta E_{bind.}^{sep.}$.
For cluster (H$_2$O)$_{5}${NH$_4$}$^+$, the five low-energy isomers are illustrated in Figure \ref{fig:nh4-4-6w}. The isomer 5-a is the most stable one, which is consistent with F. Spiegelmans result using the global Monte Carlo optimization and G. Shieldss results obtained with a mixed molecular dynamics/quantum mechanics moldel.\cite{Douady2008, Morrell2010} The energy order of 5-a to 5-e at SCC-DFTB level is consistent with that at MP2/Def2TZVP level with ZPVE correction. 5-a, 5-d and 5-e have a complete solvation shell while one dangling N-H bond is exposed in 5-b and 5-c. For the five low-energy isomers, the energy order of our results are not exactly the same with H. Changs calculation results at MP2/6-31+G(d)level with ZPVE correction.\cite{Jiang1999} In H. Changs results, 5-d is the low-energy isomer and 5-a is the second low-energy isomer. They didnt find isomers 5-b and 5-c. From the comparison, it implies the combination of SCC-DFTB and PTMD is good enough to find the low-energy isomer and the basis set can affect the energy order when using the MP2 approach.
When all the water molecules are considered as a whole part, the obtained binding energy has a deviation due to the interaction of water molecules. As listed in Table \ref{reBindE}, for isomers 5-a to 5-e, the relative binding energy $\Delta E_{bind.}^{whole}$ are -1.62, 0.72, 0.69, -1.08 and -2.08 kcal·mol\textsuperscript{-1} and $\Delta E_{bind.}^{sep.}$ are -0.56, 0.48, 0.55, -0.78 and 0.88 kcal·mol\textsuperscript{-1}, respectively. The $\Delta E_{bind.}^{whole}$ is bigger than corresponding $\Delta E_{bind.}^{sep.}$, which indicates it is better to calculate the binding energy with considering the water molecules separately. The $\Delta E_{bind.}^{sep.}$ is less than 1.00 kcal·mol\textsuperscript{-1} for the five low-energy isomers of cluster (H$_2$O)$_{5}${NH$_4$}$^+$, so the SCC-DFTB method is good enough compared to MP2/Def2TZVP with BSSE correction for cluster (H$_2$O)$_{5}${NH$_4$}$^+$.
For cluster (H$_2$O)$_{6}${NH$_4$}$^+$, no N-H bond is exposed in the five low-energy isomers displayed in Figure \ref{fig:nh4-4-6w}. 6-a is the first low-energy isomer at SCC-DFTB level, which is a symmetric double-ring species connected together by eight hydrogen bonds making it a robust structure. 6-a is also the first low-energy isomer obtained using the Monte Carlo optimizations with the intermolecular polarizable potential.\cite{Douady2008} 6-d is the first low-energy isomer at MP2/Def2TZVP level with ZPVE correction but it is only 0.22 kcal·mol\textsuperscript{-1} lower than 6-a. In Shieldss results, 6-d is also the first low-energy isomer at MP2/aug-cc-pVDZ level.\cite{Morrell2010} In H. Changs study, 6-b with a three-coordinated H2O molecule is the first low-energy isomer for cluster (H$_2$O)$_{6}${NH$_4$}$^+$ at B3LYP/6-31+G(d) level.\cite{Wang1998} 6-b is also the first low-energy isomer at B3LYP/6-31++G(d,p) level including the harmonic ZPE contribution.\cite{Douady2008} The energy of 6-b is only 0.14 kcal·mol\textsuperscript{-1} higher than that of 6-a at MP2/Def2TZVP level with ZPVE correction. The energies of 6-a, 6-b and 6-d are very close at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels, which implies it is easy to have a transformation among 6-a, 6-b and 6-d. It shows SCC-DFTB is good to find the low-energy isomers of cluster (H$_2$O)$_{6}${NH$_4$}$^+$ compared to MP2 and B3LYP methods.
As shown in Table \ref{reBindE}, for isomers 6-a to 6-e, the relative binding energy $\Delta E_{bind.}^{whole}$ are -1.71, -1.14, -2.06, -2.90 and -1.18 kcal·mol\textsuperscript{-1} and the $\Delta E_{bind.}^{sep.}$ are -0.38, -0.76, 0.27, -1.06 and -0.60 kcal·mol\textsuperscript{-1}, respectively. It indicates the binding energy are very close at SCC-DFTB and MP2/Def2TZVP with BSSE correction levels when water molecules are calculated separately. The $\Delta E_{bind.}^{whole}$ is bigger than corresponding $\Delta E_{bind.}^{sep.}$ because of the interaction of water molecules when all the water molecules are considered as a whole part.
For cluster (H$_2$O)$_{7}${NH$_4$}$^+$, the five low-energy isomers are shown in Figure \ref{fig:nh4-7-10w}. The ion core {NH$_4$}$^+$ has a complete solvation shell in isomers 7-a to 7-e. 7-a and 7-b with three three-coordinated H$_2$O molecules are the first low-energy isomers at SCC-DFTB level. In F. Spiegelmans study, 7-a is also the first low-energy isomer using the Monte Carlo optimizations with the intermolecular polarizable potential.\cite{Douady2008} 7-c is the first low-energy isomer at MP2/Def2TZVP with ZPVE correction level including three three-coordinated water molecules. 7-c is also the first low-energy isomer at B3LYP/6-31++G(d,p) level including the harmonic ZPE contribution.\cite{Douady2008} 7-e is the first low-energy isomer with three three-coordinated H$_2$O molecules at MP2/aug-cc-pVDZ level in G. Shieldss study.\cite{Morrell2010} As illustrated in Figure \ref{fig:nh4-7-10w}, the energy difference between 7-a, 7-c and 7-e at SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels are less than 0.61 kcal·mol\textsuperscript{-1} so it is possible that the first low-energy iosmer is different when different method are applied. The energy of 7-a and 7-b are the same at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels and their structures are similar, which indicates it is easy for them to transform to each other. The results for cluster (H$_2$O)$_{7}${NH$_4$}$^+$ verify the accuracy of SCC-DFTB approach.
\begin{figure}[h!]
\includegraphics[width=1.0\linewidth]{nh4-7-10w.png}
\centering
\caption{The five low-energy isomers of clusters (H$_2$O)$_{7-10}${NH$_4$}$^+$ and the associated relative energies (in kcal·mol\textsuperscript{-1}) at MP2/Def2TZVP level with (bold) and without ZPVE correction and SCC-DFTB level (italic).}
\label{fig:nh4-7-10w}
\end{figure}
As shown in Table \ref{reBindE}, for isomers 7-a to 7-e, the relative binding energy $\Delta E_{bind.}^{whole}$ are -2.95, -2.92, -2.17, -1.28 and -3.22 kcal·mol\textsuperscript{-1} and the $\Delta E_{bind.}^{sep.}$are only -0.39, -0.38, 0.09, -1.35 and -2.27 kcal·mol\textsuperscript{-1}, respectively. It indicates the binding energies of 7-a to 7-e at SCC-DFTB agree well especially for 7-a to 7-d with those at MP2/Def2TZVP with BSSE correction level when water molecules are calculated separately. When all the water molecules are regarded as a whole part, the results of SCC-DFTB are not as good as those of the MP2 with BSSE method.
For cluster (H$_2$O)$_{8}${NH$_4$}$^+$, 8-a to 8-e are the five low-energy isomers displayed in Figure \ref{fig:nh4-7-10w}. In 8-a to 8-d, the ion core {NH$_4$}$^+$ has a complete solvation shell. 8-a is the first low-energy isomer in our calculation at SCC-DFTB level. In F. Spiegelmans study, 8-b is the first low-energy isomer at B3LYP/6-31++G(d,p) level including the harmonic ZPE contribution.\cite{Douady2008} The structures of 8-a and 8-b are very similar and the energy differences are only 0.09 and 0.18 kcal·mol\textsuperscript{-1} at SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels, respectively. 8-d with seven three-coordinated H$_2$O molecules in the cube frame is the first low-energy isomer in our calculation at MP2/Def2TZVP with ZPVE correction level, which is consistent with F. Spiegelmans results obtained using Monte Carlo optimizations.\cite{Douady2008} In 8-e, {NH$_4$}$^+$ has an exposed N-H bond and it also has seven three-coordinated H$_2$O molecules in its cage frame. The energies of isomers 8-a to 8-e are very close calculated using SCC-DFTB and MP2 methods, so its possible that the energy order will change when different methods or basis sets are applied. The results certificate the SCC-DFTB is good enough to find the low-energy isomers for cluster (H$_2$O)$_{8}${NH$_4$}$^+$.
As shown in Table \ref{reBindE}, for isomers 8-a to 8-e, the relative binding energy $\Delta E_{bind.}^{whole}$are -2.20, -1.61, -3.71, -2.43 and -0.55 kcal·mol\textsuperscript{-1}, respectively and the biggest $\Delta E_{bind.}^{sep.}$ is -2.01 kcal·mol\textsuperscript{-1}. It shows the binding energies at SCC-DFTB level agree well with those at MP2/Def2TZVP with BSSE correction level when water molecules are calculated separately. From these results, when all the water molecules are considered as a whole part, the results of SCC-DFTB didnt agree well with those of the MP2 with BSSE correction method.
For cluster (H$_2$O)$_{9}${NH$_4$}$^+$, the five low-energy structures of (H$_2$O)$_{9}${NH$_4$}$^+$ are illustrated in Figure \ref{fig:nh4-7-10w}. 9-a with seven three-coordinated H$_2$O molecules in the cage frame is the first low-energy isomer at SCC-DFTB level. 9-a is also the first low-energy structure at B3LYP/6-31++G(d,p) level including the harmonic ZPE contribution in F. Spiegelmans study.\cite{Douady2008} 9-b with one N-H bond exposed in {NH$_4$}$^+$ is the second low-energy isomer whose energy is only 0.22 kcal·mol\textsuperscript{-1} higher than that of 9-a in the results of SCC-DFTB calculation. 9-b is the first low-energy isomer at MP2/Def2TZVP with ZPVE correction level in our calculation and it is also the first low-energy isomer at B3LYP/6-31++G(d,p) level in F. Spiegelmans study.\cite{Douady2008} 9-c, 9-d and 9-e have a complete solvation shell. All the water molecules are connected together in the structure of 9-c. The structures of 9-a and 9-e are very similar and their energy difference is only 0.11 kcal·mol\textsuperscript{-1} at MP2/Def2TZVP with ZPVE correction level. The energy difference of isomers 9-a to 9-e is less than 0.51 kcal·mol\textsuperscript{-1} at SCC-DFTB and less than 0.86 kcal·mol\textsuperscript{-1} at MP2/Def2TZVP with ZPVE correction, so its easy for them to transform to each other making it possible for the variation of the energy order. The results certificate the SCC-DFTB is good enough to find the low-energy isomers for cluster (H$_2$O)$_{9}${NH$_4$}$^+$.
As shown in Table \ref{reBindE}, for isomers 9-a to 9-e, the relative binding energy $\Delta E_{bind.}^{whole}$ are -2.20, -1.61, -3.71, -2.43 and -0.55 kcal·mol\textsuperscript{-1} and the relative binding energy $\Delta E_{bind.}^{sep.}$ are -1.39, -0.84, -0.85, -1.78, and -0.91 kcal·mol\textsuperscript{-1}, respectively.
It is obvious that the absolute values of $\Delta E_{bind.}^{whole}$ are bigger than the corresponding $\Delta E_{bind.}^{sep.}$. It shows the binding energies at SCC-DFTB level agree well with those at MP2/Def2TZVP with BSSE correction level when water molecules are calculated separately. According to the results, When all the water molecules are considered as a whole part, the results of SCC-DFTB didnt agree well with those of the MP2 with BSSE correction method.
For cluster (H$_2$O)$_{10}${NH$_4$}$^+$, 10-a to 10-e are the five low-energy isomers in which the ion core {NH$_4$}$^+$ has a complete solvation shell shown in Figure \ref{fig:nh4-7-10w}. 10-a with eight three-coordinated H2O molecules in its big cage structure is the first low-energy isomer calculated using the SCC-DFTB approach. 10-a is also the first low-energy structure at B3LYP/6-31++G(d,p) level including the harmonic ZPE contribution in F. Spiegelmans study.\cite{Douady2008} In 10-b and 10-e, there is a four-coordinated H$_2$O molecule in their cage structures. 10-d is the first low-energy structure in our calculation results using MP2/Def2TZVP with ZPVE correction, which is also the first low-energy isomer at B3LYP/6-31++G(d,p) level in F. Spiegelmans study.\cite{Douady2008} The energy of 10-b is only 0.17 kcal·mol\textsuperscript{-1} higher than that of 10-a at SCC-DFTB level, and it is only 0.31 kcal·mol\textsuperscript{-1} lower than that of 10-a at MP2/Def2TZVP with ZPVE correction level. The energy of isomers 10-a to 10-e are very close at both SCC-DFTB and MP2/Def2TZVP levels, which indicates the results with SCC-DFTB agree well with those using MP2/Def2TZVP method for cluster (H$_2$O)$_{10}${NH$_4$}$^+$.
As shown in Table \ref{reBindE}, for isomers 10-a to 10-e, the relative binding energies $\Delta E_{bind.}^{whole}$ and $\Delta E_{bind.}^{sep.}$ are not as small as the corresponding ones of clusters (H$_2$O)$_{1-9}${NH$_4$}$^+$, which implies the error of the relative binding energy increases with the number of water molecules in the cluster. The whole results of $\Delta E_{bind.}^{whole}$ are still bigger than those of $\Delta E_{bind.}^{sep.}$ for isomers 10-a to 10-e.
\subsubsection{Properties of (H$_2$O)$_{4-10}${NH$_3$} Clusters}
For cluster (H$_2$O)$_{4}${NH$_3$}, the five low-energy structures 4$^\prime$-a to 4$^\prime$-e are displayed in Figure \ref{fig:nh3-4-7w}. 4$^\prime$-a with three N-H bonds exposed is the first low-energy isomer at SCC-DFTB level. 4$^\prime$-b with two N-H bonds exposed is the second low-energy isomer at SCC-DFTB level but it is the first low-energy isomer at MP2/Def2TZVP with ZPVE correction level. The energy differences between 4$^\prime$-a to 4$^\prime$-b are only 0.20 and 0.07 kcal·mol\textsuperscript{-1} at SCC-DFTB and MP2/Def2TZVP with ZPVE correction level, respectively. The energy difference of isomers 4$^\prime$-a to 4$^\prime$-e is less than 0.75 kcal·mol\textsuperscript{-1} at MP2/Def2TZVP with ZPVE correction, so its possible for the variation of the energy order when different methods or basis sets are used. 4$^\prime$-d with a nearly planar pentagonal structure with nitrogen atom and the four oxygen atoms at the apexes is the first low-energy isomer at MP2/6-31+G(d,p) studied by J. Novoa et al\cite{Lee1996} 4$^\prime$-d is also the first low-energy isomer in D. Bacelos study using QCISD(T) for a single-point energy calculation based on the MP2/6-311++G(d,p) results.\cite{Bacelo2002} In addition, 4$^\prime$-a to 4$^\prime$-e are also the five low-energy isomers in D. Bacelos study even the energy order is different.\cite{Bacelo2002} The results show the SCC- DFTB is good enough to find the low-energy isomers isomers for cluster (H$_2$O)$_{4}${NH$_3$}.
\begin{figure}[h!]
\includegraphics[width=1.0\linewidth]{nh3-4-7w.png}
\centering
\caption{The five low-energy isomers of cluster (H$_2$O)$_{4-7}${NH$_3$} and the associated relative energies (in kcal·mol\textsuperscript{-1}) at MP2/Def2TZVP level with (bold) and without ZPVE correction and SCC-DFTB level (italic).}
\label{fig:nh3-4-7w}
\end{figure}
The relative binding energies of isomers 4$^\prime$-a to 4$^\prime$-e are shown in Table \ref{reBindE}. Except 4$^\prime$-d, the values of $\Delta E_{bind.}^{whole}$ for 4$^\prime$-a to 4$^\prime$-e are smaller than the corresponding values of $\Delta E_{bind.}^{sep.}$. The $\Delta E_{bind.}^{sep.}$ of 4$^\prime$-d is smaller than those of other isomers. 4$^\prime$-d has a nearly planar pentagonal structure that only contains three O-H···O hydrogen bonds among the four water molecules while other isomers contain four O-H···O hydrogen bonds among the four water molecules. So the intermolecular interaction of the four water molecules in 4$^\prime$-d is not as strong as it in other isomers, this may explain the $\Delta E_{bind.}^{sep.}$ of 4$^\prime$-d is smaller than those of other isomers. In general, both relative binding energies $\Delta E_{bind.}^{sep.}$ and $\Delta E_{bind.}^{sep.}$ are not big that indicates SCC-DFTB performs well compared to the MP2 method with BSSE correction for calculating the binding energy of cluster (H$_2$O)$_{4}${NH$_3$}.
For cluster (H$_2$O)$_{5}${NH$_3$}, 5$^\prime$-a to 5$^\prime$-e are the five low-energy isomers shown in Figure \ref{fig:nh3-4-7w}. 5$^\prime$-a with four three-coordinated water molecules is the first low-energy structure at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. 5$^\prime$-b and 5$^\prime$-c are the second and third isomers at SCC-DFTB level and they are the third and second isomers at MP2/Def2TZVP level with ZPVE. The energy difference between 5$^\prime$-b and 5$^\prime$-c is only 0.05 and 0.44 kcal·mol\textsuperscript{-1} at SCC-DFTB level and MP2/Def2TZVP with ZPVE correction level, respectively. In addition, the structures of 5$^\prime$-b and 5$^\prime$-c are very similar so it is possible for them to transform to each other. 5$^\prime$-d with two three-coordinated water molecules is the fourth low-energy structure at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. 5$^\prime$-e with four three-coordinated water molecules is the fifth low-energy structure at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. The frames of 5$^\prime$-a and 5$^\prime$-e are almost the same but the water molecule who offers the hydrogen or oxygen to form the O-H···O hydrogen bonds has a small difference. The energy of 5$^\prime$-e is 1.51 kcal·mol\textsuperscript{-1} higher than that of 5$^\prime$-a at MP2/Def2TZVP with ZPVE correction level, which implies the intermolecular connection mode has an influence on the stability of the isomers. The results show the SCC-DFTB approach performs well to find the low-energy isomers for cluster (H$_2$O)$_{5}${NH$_3$} compared with MP2/Def2TZVP with ZPVE correction method.
The relative binding energies of isomers 5$^\prime$-a to 5$^\prime$-e are shown in Table \ref{reBindE}. The values of $\Delta E_{bind.}^{whole}$ are less than 0.82 kcal·mol\textsuperscript{-1} for 5$^\prime$-a to 5$^\prime$-e. The values of $\Delta E_{bind.}^{sep.}$ are bigger than the corresponding values of $\Delta E_{bind.}^{whole}$. It indicates SCC-DFTB agrees better with MP2/Def2TZVP $\Delta E_{bind.}^{whole}$ when all the water molecules are regarded as a whole part than considering separately for calculating the binding energy of cluster (H$_2$O)$_{5}${NH$_3$}.
For cluster (H$_2$O)$_{6}${NH$_3$}, the five low-energy structures 6$^\prime$-a to 6$^\prime$-e are displayed in Figure \ref{fig:nh3-4-7w}. 6$^\prime$-a is the first low-energy structure at SCC-DFTB level. All water molecules in 6$^\prime$-a are three-coordinated. 6$^\prime$-b is the second low-energy isomer at SCC-DFTB level and its only 0.05 and 0.42 kcal·mol\textsuperscript{-1} higher than the ones of 6$^\prime$-a at SCC-DFTB level and MP2/Def2TZVP with ZPVE correction level, respectively. 6$^\prime$-c to 6$^\prime$-d are the third and fourth low-energy isomers in which the six water molecules form a triangular prism structure and there are one and two four-coordinated water molecules in 6$^\prime$-c to 6$^\prime$-d, respectively. 6$^\prime$-e is the fifth low-energy structure at SCC-DFTB level but its the first low-energy isomer at MP2/Def2TZVP with ZPVE correction level. The energy of 6$^\prime$-a to 6$^\prime$-e are very close at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels that it is difficult to keep the energy order when different methods or basis sets are applied. This also shown the SCC-DFTB method used is efficient to find the low-energy isomers of cluster (H$_2$O)$_{6}${NH$_3$}.
The relative binding energies of isomers 6$^\prime$-a to 6$^\prime$-e are listed in Table \ref{reBindE}. The smallest and the biggest values of $\Delta E_{bind.}^{whole}$ are -0.05 and -1.11 kcal·mol\textsuperscript{-1}, respectively. The smallest absolute value of $\Delta E_{bind.}^{sep.}$ is 1.96 kcal·mol\textsuperscript{-1}. The binding energies calculated with SCC-DFTB agree well with those calculated at MP2/Def2TZVP level for cluster (H$_2$O)$_{6}${NH$_3$} when all the water molecules are considered as a whole part.
For cluster (H$_2$O)$_{7}${NH$_3$}, the five low-energy isomers 7$^\prime$-a to 7$^\prime$-e are illustrated in Figure \ref{fig:nh3-4-7w}. 7$^\prime$-a with a cubic structure is the first low-lying energy structure at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. 7$^\prime$-b is the second low-energy structure at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. 7$^\prime$-b has a similar structure with 7$^\prime$-a but the NH$_3$ in it has two exposed N-H bonds. 7$^\prime$-c and 7$^\prime$-d have similar structures and they are the third and fourth low-lying energy isomers at SCC-DFTB level and their energy difference is only 0.74 kcal·mol\textsuperscript{-1}. 7$^\prime$-e with three exposed N-H bonds is the fifth low-energy isomer at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. The results of SCC-DFTB method agree well with those of MP2/Def2TZVP with ZPVE correction for the five low-energy isomers of cluster (H$_2$O)$_{7}${NH$_3$}.
The smallest and the biggest values of $\Delta E_{bind.}^{whole}$ of isomers 7$^\prime$-a to 7$^\prime$-e are -0.02 and -1.11 kcal·mol\textsuperscript{-1}, respectively and the smallest absolute value of $\Delta E_{bind.}^{sep.}$ is 2.02 kcal·mol\textsuperscript{-1} shown in Table \ref{reBindE}. The binding energies calculated with SCC-DFTB agree well with those obtained using MP2/Def2TZVP for cluster (H$_2$O)$_{7}${NH$_3$} when all the water molecules are considered as a whole part.
For cluster (H$_2$O)$_{8}${NH$_3$}, 8$^\prime$-a to 8$^\prime$-e are the five low-energy structures shown in Figure \ref{fig:nh3-8-10w}. 8$^\prime$-a in which eight water molecules constitute a cube is the first low-lying energy structure in SCC-DFTB calculation results. 8$^\prime$-b also with a water-cube structure is the second low-energy structure at SCC-DFTB level and it is the first low-energy isomer at MP2/Def2TZVP with ZPVE correction level. The energy differences between 8$^\prime$-a and 8$^\prime$-b are only 0.93 an 0.30 kcal·mol\textsuperscript{-1} at SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. From Figure \ref{fig:nh3-8-10w}, the fifth low-energy isomer 8$^\prime$-e includes less number of hydrogen bonds than other isomers and its energy has a clearly increase compared to other isomers. The results show the SCC-DFTB method performs well to obtain the low-energy isomers of cluster (H$_2$O)$_{8}${NH$_3$}.
\begin{figure}[h!]
\includegraphics[width=1.0\linewidth]{nh3-8-10w.png}
\centering
\caption{The five low-energy isomers of clusters (H$_2$O)$_{8-10}${NH$_3$} and the associated relative energies (in kcal·mol\textsuperscript{-1}) at MP2/Def2TZVP level with (bold) and without ZPVE correction and SCC-DFTB level (italic).}
\label{fig:nh3-8-10w}
\end{figure}
The smallest and the biggest values of $\Delta E_{bind.}^{whole}$ of isomers 8$^\prime$-a to 8$^\prime$-e are -0.1 and -1.28 kcal·mol\textsuperscript{-1}, respectively while the smallest absolute value of $\Delta E_{bind.}^{sep.}$ is 3.04 kcal·mol\textsuperscript{-1} shown in Table \ref{reBindE}. The binding energies calculated with SCC-DFTB agree better with those obtained at MP2/Def2TZVP level when all the water molecules are considered as a whole part in cluster (H$_2$O)$_{8}${NH$_3$} than the ones when water molecules calculated separately.
For cluster (H$_2$O)$_{9}${NH$_3$}, 9$^\prime$-a to 9$^\prime$-e are the five low-lying energy structures displayed in Figure \ref{fig:nh3-8-10w}. 9$^\prime$-a with a “chair” structure is the first low-energy structure at SCC-DFTB level. 9$^\prime$-b, 9$^\prime$-c and 9$^\prime$-d in which the nine water molecules have the similar configuration are the second, third and fourth isomers. In 9$^\prime$-b and 9$^\prime$-c, the NH$_3$ has three exposed N-H bonds and the energies of 9$^\prime$-b and 9-c are very close at both SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. The NH$_3$ has two exposed N-H bonds in 9$^\prime$-d. 9$^\prime$-e is the fifth low-energy isomer in the SCC-DFTB calculation results but it is the first low-energy isomer in the calculation results of MP2/Def2TZVP with ZPVE correction. 9$^\prime$-e has a pentagonal prism structure and all the water molecules in it are three-coordinated. The relative energy for each isomer between SCC-DFTB level and MP2/Def2TZVP with ZPVE correction level is less than 1.23 kcal·mol\textsuperscript{-1}. This shows our SCC-DFTB calculation results are consistent with the calculation results of MP2/Def2TZVP with ZPVE correction for low-energy isomers optimization of cluster (H$_2$O)$_{9}${NH$_3$}.
The relative binding energies of isomers 9$^\prime$-a to 9$^\prime$-e are shown in Table \ref{reBindE}. The absolute values of $\Delta E_{bind.}^{whole}$ are less than 1.09 kcal·mol\textsuperscript{-1} while the smallest absolute value of $\Delta E_{bind.}^{sep.}$ is 2.57 kcal·mol\textsuperscript{-1}. The binding energies calculated with SCC-DFTB agree well with those acquired at MP2/Def2TZVP level when all the water molecules are considered as a whole part for cluster (H$_2$O)$_{9}${NH$_3$}.
For cluster (H$_2$O)$_{10}${NH$_3$}, 10$^\prime$-a to 10$^\prime$-e are the five low-energy structures illustrated in Figure \ref{fig:nh3-8-10w}. The energy order for the five low-energy structures is the same at SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. 10$^\prime$-a and 10$^\prime$-b are the first and second low-energy isomer in which the ten water molecules constitute the pentagonal prism. The energy differences of 10$^\prime$-a and 10$^\prime$-b are only 0.27 and 0.58 kcal·mol\textsuperscript{-1} at SCC-DFTB and MP2/Def2TZVP with ZPVE correction levels. 10$^\prime$-c and 10$^\prime$-d are the third and fourth low-energy isomers in which eight water molecules constitute a cube and the energy differences between 10$^\prime$-c and 10$^\prime$-d are very small calculated with SCC-DFTB or MP2/Def2TZVP with ZPVE correction. 10$^\prime$-e is the fifth low-energy structure in which eight water molecules also constitute a cube but its energy is obviously higher than those of 10$^\prime$-c and 10$^\prime$-d. The calculation results of SCC-DFTB are consistent with those of MP2/Def2TZ for the optimization of the low-energy isomers of cluster (H$_2$O)$_{10}${NH$_3$}. According to the structures of the five low-energy isomers of clusters (H$_2$O)$_{1-10}${NH$_3$}, in most cases, the NH$_3$ usually contains two or three exposed N-H bonds.
The smallest and biggest values of $\Delta E_{bind.}^{whole}$ of isomers 10$^\prime$-a to 10$^\prime$-e are -0.03 and -1.10 kcal·mol\textsuperscript{-1} while the smallest absolute value of $\Delta E_{bind.}^{sep.}$ is 4.80 kcal·mol\textsuperscript{-1} shown in Table \ref{reBindE}. The values of $\Delta E_{bind.}^{whole}$ implies that SCC-DFTB agree very well with MP2/Def2TZVP for cluster (H$_2$O)$_{10}$NH$_3$ when all the water molecules are regarded as a whole part.
\subsubsection{Properties of (H$_2$O)$_{20}${NH$_4$}$^+$ Cluster}
For cluster (H$_2$O)$_{20}${NH$_4$}$^+$, the lowest-energy structure shown in Figure \ref{fig:nh3-nh4-20w} (a) was obtained with the combination of SCC-DFTB and PTMD which is consistent with that of previous study.\cite{Kazimirski2003, Douady2009, Bandow2006}
Microcanonical and canonical caloric curves were obtained using exchange Monte Carlo simulations by F. Spiegelmans group.\cite{Douady2009}
I also calculated the canonical heat capacities of cluster (H$_2$O)$_{20}${NH$_4$}$^+$ using the combination of SCC-DFTB and PTMD depicted in Figure.
\begin{figure}[h!]
\includegraphics[width=0.6\linewidth]{nh3-nh4-20w.jpeg}
\centering
\caption{The five low-energy isomers of cluster (H$_2$O)$_{20}${NH$_4$}$^{+}$ (a) and (H$_2$O)$_{20}${NH$_3$} (b) at SCC-DFTB level.}
\label{fig:nh3-nh4-20w}
\end{figure}
\subsection{Conclusions for Ammonium/Ammonia Including Water Clusters}
The low-energy isomers reported for (H$_2$O)$_{1-10, 20}${NH$_4$}$^+$ and (H$_2$O)$_{1-10, 20}$NH$_3$ clusters are obtained with a combination of
SCC-DFTB (0.12/1.16) and PTMD. Binding energies as a function of the N---O distance in (H$_2$O){NH$_4$}$^+$ and (H$_2$O)NH$_3$ demonstrate
that the improve parameters I proposed are in much better agreement with reference calculations than the original SCC-DFTB parameters.
The low-energy isomers of clusters (H$_2$O)$_{1-10}${NH$_4$}$^+$ and (H$_2$O)$_{1-10}$NH$_3$ at the SCC-DFTB (0.12/1.16) level
agree well with those at MP2/Def2TZVP level and the corresponding results in the literature. The SCC-DFTB binding energies also agree well
with those calculated with MP2/Def2TZVP method with BSSE correction. This demonstrate that SCC-DFTB (0.12/1.16) approach is good
enough to model ammonium and ammonia containing water clusters.
Among the five lowest-energy structures of (H$_2$O)$_{4}${NH$_4$}$^+$, four of them display a dangling N-H bond. Among the five lowest-energy
structures of (H$_2$O)$_{5}${NH$_4$}$^+$, only two structures display a dangling N-H bond. Among the five lowest-energy isomers of (H$_2$O)$_{6-10}${NH$_4$}$^+$,
all the structures, except 8-e, display a ion core {NH$_4$}$^+$ that has a complete solvation shell but it is not located at the center of the water cluster.
In the most stable structures of (H$_2$O)$_{20}${NH$_4$}$^+$, the ion core {NH$_4$}$^+$ has a complete solvation shell and it is in the center of the
water cluster. In contrast, in the low-energy isomers of (H$_2$O)$_{1-10}$NH$_3$ clusters, NH$_3$ is never fully solvated by the water cluster.
The present study demonstrate thet ability of SCC-DFTB to model small size ammonium and ammonia containing water clusters, which is less
expensive than \textit{ab initio} methods. It is possible for SCC-DFTB to describe the larger scaled ammonium and ammonia containing water clusters.
\section{Structural and Energetic Properties of Protonated Uracil Water Clusters} \label{structureUH}
\subsection{General introduction}
Gas phase investigations of molecules help to understand the intrinsic properties of molecules that are free from the effects of solvents.
The gas phase study needs to be extended towards more realistic biomolecular systems, to reveal how the intrinsic molecular properties are affected by the surrounding medium when the biomolecules are in a natural environment.\cite{Maclot2011, Domaracka2012, Markush2016, Castrovilli2017} The hydration study of biomolecules is of paramount importance to get insights into their behavior in aqueous medium, especially the effects on their structure, stability and dynamics.
The nucleobases in DNA and RNA play a significant role in the encoding and expression of genetic information in living systems while water is a natural medium of many reactions in living organisms. The study of the interaction between nucleobase molecules and aqueous environment has attracted a lot of interests among biologists and chemists. Exploring the clusters composed of nucleobase molecules with water is a good workbench to observe how the properties of nucleobase molecules change when going from isolated gas-phase to hydrated species.
The radiation can cause damages on RNA and DNA molecules, which is proficiently applied in radiotherapy for cancer treatment. The major drawback in radiotherapy is the unselective damage in both healthy and tumor cells, which has a big side effect. This makes it particularly important to explore the radiation fragments.
Uracil , C$_4$H$_4$N$_2$O$_2$, is one of the four nucleobases of RNA, has been paid attention concerning radiation damage. Protonated uracil UH$^+$ can be generated by radiation damages.\cite{Wincel2009}
The reasons for such degradation can be due to the interaction with slow electrons, as shown by the work of B. Boudaiffa \textit{et al.} \cite{Boudaiffa2000}
Several studies have been devoted to the effect of hydration on the electron affinity of DNA nucleobases. \cite{Smyth2011, Siefermann2011, Alizadeh2013} For instance, A. Rasmussen \textit{et al.} found that a water molecule is more likely to interact with a charged species than with a neutral one though the study of hydration effects on the lowest triplet states of cytosine, uracil, and thymine by including one or two water molecules explicitly, \cite{Rasmussen2010}
However, a lot of work is still needed to be performed to understand the role of aqueous environment on charged nucleobases of DNA and RNA.
Collision experiments is a useful tool that can be applied to understand the reactivity of molecules and provide access to structural information.\cite{Coates2018}
Fragmentation of the bare protonated U has already been performed under collision-induced dissociation (CID) with tandem mass spectrometry,\cite{Nelson1994, Sadr2014, Molina2016}
however, there are only few studies available concerning the effect of hydration on such process. Infrared photodissociation spectroscopy of singly
hydrated protonated uracil shows that the most stable tautomeric form of the neutral uracil (diketo) differs from the most stable one for bare
protonated uracil (keto-enol).\cite{Bakker2008}
%Some theoretical studies reported the structures taken by neutral microhydrated uracil, containing up to 15 water molecules.\cite{Shishkin2000, Gadre2000, Gaigeot2001, Danilov2006, Bacchus2015}
However, fragmentation studies of such species under CID conditions have not been performed. \textbf{S. Zamith and J.-M. L'Hermite conducted such CID
experiments on protonated uracil water species (H$_2$O)$_{1-15}$UH$^+$ during my thesis and I collaborated with them in order to provide a theoretical
support to their measurements}.
%The dissociation of mass selected (H$_2$O)$_{1-15}$UH$^+$ clusters is induced by collisions with target rare gas atoms (Ne, Ar) or water molecules (H$_2$O, D$_2$O) at a controlled center of mass collision energy 7.2 eV which is chosen high enough so that a large number of fragmentation channels are explored. After the collisions, the remaining charged species, which have lost one or several neutral subunits, are detected. We find essentially the same results whatever the nature of target atoms of molecules is. The resulting inter-molecular dissociation patterns of the (H$_2$O)$_{1-15}$UH$^+$ clusters show that below $n$ = 5 only water molecules are evaporated whereas, for n $\geq$ 5, a new fragmentation channel appears that corresponds to the loss of neutral uracil.
Theoretical studies have already been devoted to mixed uracil-water clusters and intended to describe the lowest energy structures. However,
only neutral species ((H$_2$O)$_{n}$U) were considered.\cite{Shishkin2000, Gadre2000, Van2001diffu, Gaigeot2001, Danilov2006, Bacchus2015}
Those studies showed that for sizes up to with $n$ = 3, the water molecules arrange in monomers or dimers in the plane of the uracil molecule
\cite{Gadre2000, Van2001diffu, Gaigeot2001, Danilov2006, Bacchus2015} with no trimer formation. But for $n$ \textgreater~3, very different structures
were predicted depending on the considered study. For instance, M. Ghomi predicted that for $n$ = 7,\cite{Gaigeot2001} water molecules arrange
in dimers and trimers in the plane of the uracil molecule, whereas for $n$ = 11, water molecules form locked chains.\cite{Shishkin2000} 3D configurations were also proposed. For instance, all water molecules lie above the uracil plane for $n$ = 4, 5 reported by F. Calvo \textit{et al.}.\cite{Bacchus2015} Similarly, for $n$ = 11, V. Danilov \textit{et al.} also obtained a structure that consists of a water cluster above the uracil molecule.\cite{Danilov2006} Such structures are predicted to start with 4 water molecules
reported by F. Calvo and collaborator \cite{Bacchus2015} or with 6 water molecules (though 5 have not been calculated) reported by S. Gadre \textit{et al.}.\cite{Gadre2000}
Those studies may suggest that for few water molecules (up to two), the proton should be located on the uracil molecule, whereas when a large number
of water molecules surround the uracil, the charge is expected to be located on the water molecules. Of course, the excess proton is expected to strongly
influence the structure of the lowest energy isomers of each species, as observed for pure water clusters, so the size at which the proton is transferred
from uracil to water cannot be deduced from the aforementioned studies. Moreover, all those theoretical studies do not lead to the same low-energy
structures as highlighted by V. Danilov and F. Calvo.\cite{Danilov2006, Bacchus2015} Consequently, although it is instructive from a qualitative point
of view, the analysis of the experimental data by S. Zamith and J.-M. L'Hermite cannot be based on those studies. I have therefore undertaken a
theoretical simulation of hydrated protonated uracil clusters (H$_2$O)$_{1-7, 11, 12}$UH$^+$ to determine their lowest-energy structures to
complete the experiments by S. Zamith and J.-M. L'Hermite at the \textit{Laboratoire Collisions Agr\'egats R\'eactivit\'e }(LCAR). This work has
been published in 2019 in the \textit{The Journal of Chemical Physics}.\cite{Braud2019}
\subsection{Results and Discussion}
In the following section, section~\ref{exp_ur}, I present in details the results obtained from the CID experiments of S. Zamith and J.-M. L'Hermite
and the main concepts used to interpret the data. The following section, section~\ref{calcul_ur}, is devoted to the theoretical determination of
the low-energy isomers of the (H$_2$O)$_{1-7,11,12}$UH$^+$ clusters. A more detailed presentation of CID experiments is also provided in
section~\ref{exp_cid} of chapter~4, where these details are important to explicitly model CID experiments.
\subsubsection{Experimental Results} \label{exp_ur}
\textbf{Time-of-flight of mass spectrum.}
A typical fragmentation mass spectrum obtained by colliding (H$_2$O)$_{7}$UH$^+$ with neon at a center of mass collision energy of 7.2 eV is shown in Figure \ref{mass7w}. The more intense peak on the right comes from the parent cluster (H$_2$O)$_{7}$UH$^+$, the next 7 peaks at the left of the parent peak correspond to the loss of 1-7 water molecules of parent cluster, and the next 5 peaks to the left results from the evaporation of the uracil molecule and several water molecules from parent cluster. This mass spectrum is obtained at the highest pressure explored in the present experiments. This is still true for the largest size investigated here, namely, (H$_2$O)$_{15}$UH$^+$. From the result of the fragmentation mass spectrum displayed in Figure \ref{mass7w}, it indicates multiple collisions are possible, which allows the evaporation of all water molecules. Moreover, the intensity of evaporation of water molecules is bigger than the one of evaporation of U.
In the study, mainly focus on two specific channels. Channel 1 corresponds to the loss of only neutral water molecules, whereas channel 2 corresponds to the loss of neutral uracil and one or several water molecules,
\begin{align}
\mathrm{Channel~1, ~(H_2O)_nUH^+} & \rightarrow ~ \mathrm{(H_2O)_{n-x}UH^+ + xH_2O } \\
\mathrm{Channel~2, ~(H_2O)_nUH^+} & \rightarrow ~ \mathrm{ (H_2O)_{n-x}H^+ + xH_2O + U}
\end{align}
\figuremacro{mass7w}{Time-of-flight of mass spectrum obtained by colliding (H$_2$O)$_{7}$UH$^+$ with Ne
at 7.2 eV center of mass collision energy (93.5 eV in the laboratory frame).}
\textbf{Fragmentation cross section.}
The total fragmentation cross sections of clusters
(H$_2$O)$_{n-1}$UH$^+$, pure water clusters \newline (H$_2$O)$_{2-6}$H$^+$,\cite{Dalleska1993} and deuterated water clusters (D$_2$O)$_{5, 10}$H$^+$ \cite{Zamith2012} are plotted in Figure \ref{fragcrosssec} as a function of the cluster size n. Here n stands for the total number of molecules when the cluster includes uracil molecule. Different target atoms and molecules were used in these experiments: Water molecules or neon atoms in our experiments, xenon atoms in Dalleskas experiments. These experimental data are compared to the geometrical (\textit{i.e.}, hard sphere) cross sections given by:
\begin{align}
\label{cross-section-geo}
\sigma_{geo} = \pi \left(\left[n_w\times r_w^3 + n_Ur_U^3\right]^{1/3} + r_T \right)^2
\end{align}
where $n_w$ is the number of water molecules, and $n_U$ is the number of uracil molecules ($n_U$ = 0 or 1 in the present study). $r_w$, $r_U$, and $r_T$ refer to the molecular radii of water, uracil, and the target atom or molecule, respectively. The molecular radii are deduced from macroscopic densities that gives $r_U$ = 3.2 \AA \cite{Myers2007} and $r_w$ = 1.98 \AA. The radii of rare gas target atoms are taken as their Van der Waals radii $r_{Ne}$ = 1.54 \AA and $r_{Xe}$ = 2.16 \AA.
The main differences between the curves in Figure \ref{fragcrosssec} can be rationalized as follows: The larger the size of the target atom (or molecule) is, the bigger the fragmentation cross section will be. The experimental fragmentation cross sections of clusters (H$_2$O)$_{n-1}$H$^+$ colliding with water molecules are larger than the values obtained for collisions with Ne atoms. In the same vein, for a given number of molecules in the cluster, the cross section is larger for clusters containing uracil. The overall trend of all curves in Figure \ref{fragcrosssec} is the same: The
fragmentation cross sections increase with the size and seem to tend toward the geometrical one.
The cross sections measured for clusters containing uracil colliding with water molecules (black squares) are of the same magnitude as the ones previously obtained for deuterated pure water clusters (green full circles) at a similar collision energy.\cite{Zamith2012} For clusters containing uracil, fragmentation cross sections are systematically larger than the one for pure water clusters by an amount of the same magnitude as the one predicted by the geometrical cross sections. For instance, the difference between red squares and blue stars, and the difference between red full line and blue dashed line has the same magnitude.
The fragmentation cross sections obtained by N. Dalleska and coworkers \cite{Dalleska1993} for protonated water clusters are within our error bars for $n$ = 5, 6 and about a factor of 2 lower for $n$ = 3, 4. However their cross section is notably lower for (H$_2$O)$_2$H$^+$ as compared to our measurement for (H$_2$O)UH$^+$. This difference may be explained by the fact that UH$^+$ forms a weaker bond with water than H$_2$OH$^+$ does. Indeed the dissociation energy D[H$_2$OH$^+$H$_2$O] is 1.35 eV \cite{Dalleska1993, Hansen2009} whereas the value for D[UH$^+$H$_2$O] is estimated between 0.54 \cite{Wincel2009} and 0.73 eV. \cite{Bakker2008} The same behavior is observed for $n$ = 3, and the dissociation energy D[(H$_2$O)$_2$H$^+$H$_2$O] = 0.86 eV \cite{Dalleska1993, Hansen2009} is greater than the dissociation energy D[U(H$_2$O)H$^+$H$_2$O] = 0.49 eV.\cite{Wincel2009} Hence the dissociation of water molecules is more favored in the protonated uracil cluster than in the pure water clusters.
\figuremacro{fragcrosssec}{Fragmentation cross sections of clusters (H$_2$O)$_{n-1}$UH$^+$ at a collision energy of 7.2 eV plotted as a function of the total number n of molecules in the clusters. The experimental results and geometrical cross sections are shown for collision with H$_2$O and Ne. The results from N. Dalleska et al.\cite{Dalleska1993} using Xe as target atoms on pure protonated water clusters (H$_2$O)$_{2-6}$H$^+$ and from S. Zamith \textit{et al.} \cite{Zamith2012} using water as target molecules on deuterated water clusters (D$_2$O)$_{5,10}$H$^+$ are also shown. The geometrical collision cross sections of water clusters in collision with Xe atoms and water molecules are also plotted. Error bars represent one standard deviation.}
\textbf{Intermolecular fragmentation.}
Figure \ref{Uloss} displays the percentage of the fragments that have lost a neutral uracil molecule over all the fragments, plotted as a function of the number of water molecules in the parent cluster (H$_2$O)$_{n}$UH$^+$. It shows that for the cluster (H$_2$O)$_{n}$UH$^+$ with a small number of water molecules, almost no neutral uracil is evaporated. From $n$ = 5 and more clearly from $n$ = 6, the loss of neutral uracil molecule increases up to about 20\% for (H$_2$O)$_{9}$UH$^+$.
\figuremacro{Uloss}{Proportion of neutral uracil molecule loss plotted as a function of the number of water molecules n in the parent cluster (H$_2$O)$_{n}$UH$^+$. Results obtained for collisions with Ne atoms at 7.2 eV center of mass collision energy.}
The fragmentation can arise from two distinct mechanisms (direct and statistical fragmentation processes) depending on the life time of the collision complex. On the one hand, if the fragmentation occurs in a very short time after collision, the dissociation is impulsive (direct). In this case, I thus assume that the nature of the collision products is partly determined by the nature of the lowest-energy isomers of parent clusters and especially by the location of the excess proton in the structure. In other words, the lowest-energy isomer of the parent cluster obviously plays a major role in determining the fragmentation channels. On the other hand, in the case of long-lived collision complexes, collision energy is transferred to the parent cluster and is redistributed among all degrees of freedom. This is a slow process, and the structures involved during the fragmentation are no longer the lowest-energy isomer, \textit{i.e.}, the structure of the cluster can undergo structural reorganizations before evaporation. Furthermore, the excess proton can also diffuse in the structure and for instance, recombine with the uracil. Then the role of the initial structure of the parent clusters is strongly reduced in determining the fragmentation channels.
In Figure \ref{Uloss}, I focus on the loss of the neutral uracil molecule in the detected fragments since it indicates where the proton lies after collision, namely, on the uracil or on a water cluster. A transition in the nature of fragmentation product is clearly seen from $n$ = 5-6. To account for this transition, the evaporation originates from a direct fragmentation process is considered. A short discussion about the implications of possible structural rearrangement prior to dissociation, which occurs in a statistical process, will be provided in section~\ref{calcul_ur}.
The relative proton affinities of each component of the
mixed clusters gives a first estimate of which molecule, uracil or water, is more likely to carry the positive charge prior to collisions. Experimentally, the gas phase proton affinity of uracil is bracketed to 9 $\pm$ 0.12 eV.\cite{Kurinovich2002} For the proton affinity of water molecule, an experimental value is reported at 7.31 eV \cite{Magnera1991} and a theoretical one at 7.5 eV.\cite{Cheng1998} In the work of H. Cheng, it shows that the proton affinity of water clusters increases with their size.\cite{Cheng1998} The proton affinities extracted from the different studies for the uracil molecule and for water clusters as a function of the number of water molecules are displayed in Figure~\ref{protonAffinity}.
\figuremacro{protonAffinity}{The proton affinities of water clusters as a function of the number of water
molecules n, which are taken from the work of T. Magnera (black circles) \cite{Magnera1991} and from the work of Cheng (blue squares).\cite{Cheng1998} The value of the proton affinity of uracil (red dotted dashed line) is also plotted.\cite{Kurinovich2002}}
It clearly shows that the proton affinity of uracil, PA[U], is larger than the one of water monomer PA[H$_2$O]. Thus, for the mono-hydrated uracil, from the energetic point of view, the proton is on the uracil molecule and the only observed fragments are indeed protonated uracil molecules. Moreover, an experimental work \cite{Bakker2008} confirms that there is no proton transfer from the uracil to the water molecule in mono-hydrated clusters. Proton affinity of the uracil molecule is also larger than that of the water dimer, or even the trimer: PA[U] $>$ PA[(H$_2O$)$_n$], $n$ = 2 or 3 depending on the considered data for water. This is still consistent with our experimental observation of no neutral uracil molecule loss for $n$ = 2 and 3. However from the PA values, one would predict that the appearance of neutral uracil should occur for n $\approx$ 3-4. For instance, for $n$ = 4, assuming a statistical fragmentation for which the energies of final products are expected to be of relevance, the channel U + (H$_2$O)$_4$H$^+$ is energetically favorable. If one now assumes a direct dissociation, where the parent protonation state remains unchanged, one also expects that neutral uracil evaporates. However, experimentally, for $n$ = 4, no neutral uracil evaporation is observed. The loss of neutral uracil starts at $n$ = 5 and becomes significant only at $n$ = 6.
This analysis based on PA is however quite crude. Indeed, it assumes that the protonated uracil cluster would be composed of a uracil molecule attached to an intact water cluster. However, one expects that the hydration of uracil may be more complicated than this simple picture. Therefore, the uracil hydration is explored theoretically in the next Section, Section~\ref{calcul_ur}, in order to determine the proton location more realistically.
\subsubsection{Calculated Structures of Protonated Uracil Water Clusters} \label{calcul_ur}
As discussed in section~\ref{sec:ammoniumwater}, I have proposed a modified set of NH parameters to describe sp$^3$ nitrogen atoms. For,
sp$^2$ nitrogen atoms there is no need to modified the integral parameters as SCC-DFTB describe them rather correctly. Consequently, only the
$D_{NH}$ parameter needs to be defined for the present calculations. Table~\ref{tab:DNH} present the binding energy of the two
(H$_2$O)U isomers represented in Figure~\ref{uracil_i} at MP2/Def2TZVP and SCC-DFTB levels of theory. Both $D_\textrm{NH}$ = 0.12 and
$D_\textrm{NH}$ = 0.14 lead to consistent binding energies. So, to be consistent with the work performed in the previous section, we
have used $D_\textrm{NH}$ = 0.12 in the following.
\begin{figure}[h!]
\includegraphics[width=0.5\linewidth]{a-b.pdf}
\centering
\caption{Structure of two (H$_2$O)U isomers used for binding energy calculations.}
\label{uracil_i}
\end{figure}
\begin{table}
%\footnotesize
\begin{center}
%\begin{adjustbox}{width=1\textwidth}
% \small
\caption{Binding energy of two (H$_2$O)U isomers at MP2/Def2TZVP and SCC-DFTB levels of theory.}
\label{tab:DNH}
\resizebox{1.0\textwidth}{!}{
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
\textbf{isomer} & \boldm{$E_{bind_{MP2}}$} & \boldm{$E_{bind_{DFTB}}$}
& \boldm{$E_{Re}$} &\boldm{$E_{bind_{DFTB}}$}
& \boldm{$E_{Re}$}
&\boldm{$E_{bind_{DFTB}}$}
& \boldm{$E_{Re}$} \\
& & \boldm{$D_{{NH}_{0.0}}$} & \boldm{$D_{{NH}_{0.0}}$} & \boldm{$D_{{NH}_{0.12}}$} & \boldm{$D_{{NH}_{0.12}}$} & \boldm{$D_{{NH}_{0.14}}$} &\boldm{$D_{{NH}_{0.14}}$} \\
\hline
\bf{a} & -8.3 & -8.6 & -0.3 & -9.8 &-1.5 & -10.0 &-1.7 \\
\hline
\bf{b} & -6.9 & -6.6 & 0.3 & -6.9 & 0.0 & -6.9 & 0.0 \\
\hline
\end{tabular}}
% \end{adjustbox}
\end{center}
\end{table}
%\figuremacro{a-b}{The structure of isomer a and b of cluster (H$_2$O)U.}
The lowest-energy isomers determined theoretically for
hydrated uracil protonated clusters (H$_2$O)$_{1-7, 11, 12}$UH$^+$ are shown in Figures \ref{1a-f}-\ref{12a-f}. In the experiments, clusters are produced at a temperature of about 25 K, so only a very few isomers are likely to be populated. Indeed, the clusters are produced in the canonical ensemble at the temperature $T_\mathrm c \approx$ 25 K, so only isomers for which the Boltzmann factor exp(-$\Delta E k_\mathrm{B} T_\mathrm{c}$) is larger than 10$^{-7}$ are considered here. In this formula, $\Delta E$ represents the relative energy of a considered isomer with respect to the lowest-energy one. Thus for each isomer, only the six lowest-energy structures of U(H$_2$O)$_{1-7, 11, 12}$UH$^+$ obtained from the PES exploration will be discussed.
Figure \ref{1a-f} displays the six lowest-energy isomers obtained for (H$_2$O)UH$^+$. Two (1a and 1b) of them contain the u138-like isomer of U (each one with a different orientation of the hydroxyl hydrogen). Three of them (1c, 1d, and 1e) contain the u178 isomer and 1f contains the u137\cite{Wolken2000} isomer with a reverse orientation of the hydroxyl hydrogen. From those isomers, different sites are possible for the water molecule attachment which leads to variety of isomers even for such small size system. To the best of our knowledge, (H$_2$O)UH$^+$ is the most studied protonated uracil water cluster and our results are consistent with previous
published studies. Indeed, S. Pedersen and co-workers \cite{Pedersen2014} conducted ultraviolet action spectroscopy on (H$_2$O)UH$^+$ and discussed their measurements in the light of theoretical calculations performed on two isomers: ur138w8 (1a in the present study) and ur178w7 (1c).\cite{Pedersen2014} Their energy ordering at 0 K is the same whatever the computational method they used: B3LYP/6-311++G(3df,2p), M06-2X/6-311++G(3df,2p), MP2/6-311++G(3df,2p), CCSD(T)/6-311++G(3df,2p), and CCSD(T)/augcc-pVTZ and is similar to what I obtained. Similarly, J. Bakker and co-workers\cite{Bakker2008} considered three isomers: U(DK)H$^+_W$ (1a), U(KE)H$^+_{Wa}$ (1c), and U(KE)H$^+_{Wb}$ (1e) at the B3LYP/6-311++G(3df,2p) level of theory and obtained the same energy ordering as I did. Our methodology has thus allowed us to retrieve those isomers and to locate two new low-energy structures (1b and 1d). 1f is too high in energy to be considered in low-temperature experiments that are in the same range of relative energies but have never been discussed. To ensure that they are not artificially favored in our computational method, calculations were also performed at the B3LYP/6-311++G(3df,2p) level of theory. The results are presented in Figure \ref{1a-f-b3lyp}, which are consistent with the MP2/Def2TZVP ones. This makes us confident in the ability of the present methodology to locate meaningful low energy structures. Importantly, no isomer with the proton on the water molecule was obtained, neither at the DFTB or MP2 levels.
\figuremacrob{1a-f}{Lowest-energy structures of (H$_2$O)UH$^+$ obtained at the MP2/Def2TZVP level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. Important hydrogen-bond distances are indicated in bold and are given in \AA.}
\figuremacrob{1a-f-b3lyp}{Lowest-energy structures of (H$_2$O)UH$^+$ obtained at the B3LYP/6-311++G(3df,2p) level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. The corresponding values with ZPVE corrections are provided in brackets. Important hydrogen-bond distances are indicated in bold and are given in \AA.}
Figures \ref{2a-f} and \ref{3a-f} display the six lowest-energy isomers obtained for (H$_2$O)$_2$UH$^+$ and (H$_2$O)$_3$UH$^+$, respectively. For (H$_2$O)$_2$UH$^+$, the lowest energy structure, 2a contains the u138 isomer of uracil. 2b, 2d, and 2e contain u178 and 2c contains u138 with reverse orientation of the hydroxyl hydrogen. 2f contains u178 with reverse orientation of the hydroxyl hydrogen. This demonstrates that, similarly to (H$_2$O)UH$^+$, a diversity of uracil isomers are present in the low-energy structures of (H$_2$O)$_2$UH$^+$ which makes an exhaustive exploration of its PES more difficult. The same behavior is observed for (H$_2$O)$_3$UH$^+$. The configuration of u138 does not allow for the formation of a water dimer which leads to two unbound water molecules in 2a. By contrast, a water-water hydrogen bond is observed for 2b and 2c. The existence of a water dimer was not encountered in the low-energy isomers of the unprotonated (H$_2$O)$_2$U species due to the absence of the hydroxyl group on U. It is worth pointing out that 2a, 2b, 2c, and 2d are very close in energy which makes their exact energy ordering difficult to determine. However, no isomer displaying an unprotonated uracil in the low-energy isomers of (H$_2$O)$_2$UH$^+$ was located. The lowest-energy structure of (H$_2$O)$_3$UH$^+$, 3a, is characterized by two water-water hydrogen bond that forms a linear water trimer. Higher energy isomers display only one (3b, 3d, and 3e) or zero (3c and 3f) water-water bond (see Figure \ref{3a-f}). Similarly to (H$_2$O)$_2$UH$^+$, no isomer displaying an unprotonated uracil was located for (H$_2$O)$_3$UH$^+$.
\figuremacrob{2a-f}{Lowest-energy structures of (H$_2$O)$_2$UH$^+$ obtained at the MP2/Def2TZVP level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. Important hydrogen-bond distances are indicated in bold and are given in \AA.}
\figuremacrob{3a-f}{(H$_2$O)$_3$UH$^+$ lowest-energy structures obtained at the MP2/Def2TZVP level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. Important hydrogen-bond distances are indicated in bold and are given in \AA.}
The six lowest-energy isomers obtained for (H$_2$O)$_4$UH$^+$ and (H$_2$O)$_5$UH$^+$ are displayed in Figures \ref{4a-f} and \ref{5a-f}, which constitute a transition in the behavior of the proton. Indeed, in (H$_2$O)$_4$UH$^+$, two kind of low-lying energy structures appear: (i) structures composed of UH$^+$, one water trimer, and one isolated water molecule (4b, 4d, 4e, and 4f); (ii) structures composed of U and a protonated water tetramer (4a and 4c). In the latter case, the hydronium ion is always bounded to an uracil oxygen atom. The UH$_2$OH$^+$ bond is always rather strong as compared to UH$_2$O bonds as highlighted by the corresponding short oxygen-hydrogen distance. Furthermore, speaking of distances, the difference between the UH$_2$OH$^+$ and UH$^+$H$_2$O forms is rather fuzzy and might be sensitive to computational parameters and also to quantum fluctuations of the hydrogen. This suggests that collision with (H$_2$O)$_4$UH$^+$ is more likely to induce evaporation of H$_2$O rather than H$_2$OH$^+$ or a protonated water cluster. The picture is significantly different in (H$_2$O)$_5$UH$^+$ where the lowest-energy structure displays a hydronium ion separated by one water molecule from U. Such structures do not appear in (H$_2$O)$_4$UH$^+$ due to the limited number of water molecules available to separate H$_2$OH$^+$ from U. Such separation suggests that, if considering a direct dissociation process, evaporation of neutral uracil can now occurs in agreement with the experimental observations (see discussion above). One can see that 5b, which is only 0.3 kcal.mol$^{-1}$ higher in energy than 5a, still displays a UH$_2$OH$^+$ link. This is in line with the low amount of neutral uracil that is evaporated in the experiment (see Figure \ref{Uloss}).
\figuremacrob{4a-f}{Lowest-energy structures of (H$_2$O)$_4$UH$^+$ obtained at the MP2/Def2TZVP level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. Important hydrogen-bond distances are indicated in bold and are given in \AA.}
\figuremacrob{5a-f}{Lowest-energy structures of (H$_2$O)$_5$UH$^+$ obtained at the MP2/Def2TZVP level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. Important hydrogen-bond distances are indicated in bold and are given in \AA.}
Figures \ref{6a-f} and \ref{7a-f} display the six lowest-energy isomers obtained for (H$_2$O)$_6$UH$^+$ and (H$_2$O)$_7$UH$^+$. Similarly to (H$_2$O)$_5$UH$^+$, the first lowest-energy structure, 6a and 7a, for both species (H$_2$O)$_6$UH$^+$ and (H$_2$O)$_7$UH$^+$ have the excess proton on a water molecule that is separated by one water molecule from the uracil. This appears to be common to the clusters with at least 5 water molecules. This is also observed for higher-energy isomers (6c, 6d, 7c, 7e, and 7f). Other characteristics of the proton are also observed: proton in a similar Zundel form \cite{Zundel1968} bounded to the uracil (6b, 6e, and 7d) or H$_2$OH$^+$ still bounded to uracil (6f and 7b).
Finally, due to the neutral uracil loss proportion starts to decrease from $n$=9 (see Figure \ref{Uloss}), which attracted us to perform the optimization of big cluster (H$_2$O)$_{11, 12}$UH$^+$ as examples to explore why it has this change. The six low-lying energy isomers obtained for cluster (H$_2$O)$_{11, 12}$UH$^+$ are shown in Figures \ref{11a-f} and \ref{12a-f}.
In all isomers (11a to 11f) of cluster (H$_2$O)$_{11}$UH$^+$, the excess is on the water cluster and was separated by water molecule from uracil.
For 12a, 12b, 12c, and 12d, it is obvious that the excess proton is not directly bounded to the uracil. The uracil in 12a and 12d belongs to the di-keto form (there is a hydrogen atom on each nitrogen of uracil), and the excess proton was separated by one water molecule from uracil, additionally, the uracil is surrounded by the water cluster, all of these may lead the excess proton to go to the near oxygen atom of uracil. For 12b, the excess proton is on the water cluster and is very far from the uracil. For 12c, the excess proton was separately by one water molecule from uracil. For isomers 12e and 12f, the excess proton is between the uracil and a water molecule. The uracil is surrounded by the water cluster in 12e but it is not in 12f. Of course, for (H$_2$O)$_{11}$UH$^+$, (H$_2$O)$_{12}$UH$^+$, (H$_2$O)$_6$UH$^+$ and (H$_2$O)$_7$UH$^+$ and also (H$_2$O)$_4$UH$^+$ and (H$_2$O)$_5$UH$^+$, the amount of low-energy isomers is expected to be very large and do not intended to find them all. Furthermore, due to the limited number of MP2 geometry optimization I performed, there are few chances that I located the global energy minima for (H$_2$O)$_6$UH$^+$, (H$_2$O)$_7$UH$^+$, (H$_2$O)$_{11}$UH$^+$ and (H$_2$O)$_{12}$UH$^+$. However, the general picture I am able to draw from the present discussed structures fully supports the experimental results: from (H$_2$O)$_5$UH$^+$, it exists low-energy structures populated at very low temperature in which the excess proton is not directly bound to the uracil molecule. Upon fragmentation, this allows the proton to remain bounded to the water molecules.
\figuremacrob{6a-f}{Lowest-energy structures of (H$_2$O)$_6$UH$^+$ obtained at the MP2/Def2TZVP level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. Important hydrogen-bond distances are indicated in bold and are given in \AA.}
\figuremacrob{7a-f}{Lowest-energy structures of (H$_2$O)$_7$UH$^+$ obtained at the MP2/Def2TZVP level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. Important hydrogen-bond distances are indicated in bold and are given in \AA.}
\figuremacrob{11a-f}{Lowest-energy structures of (H$_2$O)$_{11}$UH$^+$ obtained at the MP2/Def2TZVP level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. Important hydrogen-bond distances are indicated in bold and are given in \AA.} \\
\figuremacrob{12a-f}{Lowest-energy structures of (H$_2$O)$_{12}$UH$^+$ obtained at the MP2/Def2TZVP level of theory. Relative ($E_\textrm{rel}$) and binding energies ($E_\textrm{bind}$) are given in kcal.mol$^{-1}$. Important hydrogen-bond distances are indicated in bold and are given in \AA.}
All the aforementioned low-lying energy structures are relevant to describe the \newline (H$_2$O)$_{1-7, 11, 12}$UH$^+$ species at low temperature and to understand the relation between the parent cluster size and the amount of evaporated neutral uracil in the case of direct dissociation. However, as already stated, one has to keep in mind that upon collision statistical dissociation can also occur. In that case, structural rearrangements are expected to occur which are important to understand each individual mass spectra of the (H$_2$O)$_{1-15}$UH$^+$ clusters and the origin of each collision product. For instance, the fragment UH$^+$ is detected for all cluster sizes in experiment. This means that for the largest sizes, for which I have shown from the calculation that the proton is located away from the uracil, proton transfer does occur prior to dissociation. One possible scenario is that after collision, water molecules sequentially evaporates. When the number of water molecules is small enough, the proton affinity of uracil gets larger than the one of the remaining attached water cluster. Proton transfer is then likely and therefore protonated uracil can be obtained at the end.
If one turns to the neutral uracil evaporation channel, it appears that the smaller clusters H$_2$OH$^+$ and (H$_2$O)$_2$H$^+$ are not present in the time-of-flight mass spectra. This absence might have two origins. First, the dissociation energies of the protonated water monomers and dimers are substantially higher than larger sizes, and they are therefore less prone to evaporation. Second, as already mentioned, for such small sizes, the proton affinity of uracil gets larger than the one of the water dimer or trimer and proton transfer to the uracil is likely to occur.
In order to confirm the above scenarios, simulations
and/or evaporation rate calculation would have to be conducted to describe the fragmentation channels in details. MD simulations of protonated uracil have already been performed by R. Spezia and co-workers to understand the collision-induced dissociation.\cite{Molina2015, Molina2016} Although, in the present case, the initial position of the excess proton appears as a key parameter to explain the evaporation of neutral uracil, such MD simulations could be additionally conducted to provide a clearer picture on the various evaporation pathways, which will be shown in section \ref{sec:collisionwUH}.
\subsection{Conclusions on (H$_2$O)$_{n}$UH$^+$ clusters}
The work in this section presents the collision-induced dissociation of hydrated protonated uracil (H$_2$O)$_{1-15}$UH$^+$ clusters and their experimental
absolute fragmentation cross sections. The experiments demonstrate that the evaporation channels evolve with size: Below $n$ = 5, the observed charged fragments
always contain the uracil molecule, whereas from $n$ = 5, the loss of a neutral uracil molecule becomes significant. To understand this transition, I conducted an
exhaustive exploration of the potential energy surface of (H$_2$O)$_{1-7, 11, 12}$UH$^+$ clusters combining a rough exploration at the SCC-DFTB level with
fine geometry optimizations at the MP2 level of theory. Those calculations show that below $n$ = 5, the excess proton is either on the uracil or on a water molecule
directly bound to uracil, \textit{i.e.}, forming a strongly bound UH$_2$OH$^+$ complex. From $n$ = 5 and above, clusters contain enough water molecules to allow
for a net separation between uracil and the excess proton: The latter is often found bound to a water molecule which is separated from uracil by at least one other
water molecule. Upon direct dissociation, the excess proton and the uracil can thus belong to different fragments. This study demonstrates that combination of
collision-induced dissociation experiments and theoretical calculation allows to probe the solvation and protonation properties of organic molecules such as nucleobases.
This is a step toward a better understanding of the role of water in the chemistry of in vivo DNA and RNA bases. However, the knowledge of the lowest-energy isomers
of the species involved in CID experiments is not enough to understand all the collision process. To get a deeper understanding of the collision mechanism, an explicit
modelling of the collision is needed. This question is addressed in the next chapter of this thesis.
% ---------------------------------------------------------------------------