As everything in life, this work is the result of a series of direct and indirect events, which have involved a certain amount of people, each one having a different weight in this thesis. If it were possible to perceive all the detailed interactions that have led to the completion of this Ph. D., from the foremost to the impercectible ones, the list of the involved people would be boring and surely useless. In despite of this, I am grateful also to all those do not know about my existence and indirectly intervened in my life up to this point.

Of course, I want to thank two cities: Messina and Paris. In these two places I have spent the most part of the last 40 months and, there, I have learned part of the modern way to do science. Indeed, in these two cities there are some institutions in which day by day I have had the opportunity to discover my brain borders and thus to expand them. Also for this reason I thank the Università degli Studi di Messina, the Université Pierre et Marie Curie -Sorbonne Université -in Paris, and the Institute for Chemical-Physical Processes of the Italian Research Council in Messina. I am grateful to these institutions mainly because of the people that (hardly) work there. The people to whom I am most tight from the scientific and human perspectives because of the importance that they had in my last-three-years growth are surely Prof. A. Marco Saitta, Dr. Franz Saija, Prof. Paolo V. Giaquinta, Prof. Fabio Pietrucci, and Prof. Salvatore Savasta. Each of them played a central and irreplaceable role in helping me to tend to be a good scientist. Prof. Giaquinta has taught me not only how the same sentence can be posed in a colorful variety of ways, but he also has shown me the best one. Prof. Savasta has displayed me the loftiness of giving without requiring something, peculiarity that is more unique than rare, in particular in Messina.

I want to spend two (maybe useless) words about Prof. Saitta, Dr. Saija, and Prof. Pietrucci.

I consider the former two a sort of older brothers. I have shared with them almost every experience that characterized my life in the last four years. Marco and Franz have deeply influenced my thought. They have shown me the way of moving, talking, thinking, interacting with science and scientists without any kind of filter. I have learned from them something that none well-payed instruction can give you. I will be in debt with them forever. Prof. Pietrucci, since his arrival in Paris in September 2014, has been for me extremely important. Fabio has had the capability to teach to a perfect ignorant (me) some methods that I never thought to learn and apply. His professionalism and precision have impressed me in a non-reversible way.

From the technical point of view, this work would be inexistent without the French National Supercomputing Facility GENCI-IDRIS which gave to me enough (a lot) of time computing.

It is straightforward to thank Sara Laporte and Adrien Mafety that have been my first friends in Paris. With them I have shared not only a good amount of reciprocal problems that a Ph. D. student have to tackle, but also a good amount of fun. A special thank goes also to Guilherme Ribeiro with which I felt often in resonance from the human perspective. In the same way, a great thank goes to all my colleagues in Messina.

Of course, I cannot forget my family assistance, from my grandparents to my brothers. My dad and my mother characterized my propensity to the scientific method as a normal way to think whereas, at the same time and together with all the other familiar figures, they have taught me the beauty of life without any lesson. My gratefulness goes to my parents also because of their "materialist" contribution. Finally, a special remark goes to Katia that constantly has assisted me with her precious presence, also "at distance", during every instant of this work, giving me a not negligible amount of vital energy.

Contents

References Introduction

The simplest mechanism by which a rearrangement of covalent bonds occurs is the ionization of a molecule through the loss/gain of a proton H + . This phenomenon is of paramount importance since it represents the foundation of the concept of pH, of the functioning of many ecological devices which exploit proton transfer (PT) (i.e., Nafion membranes, Direct Methanol Fuel Cells, etc.), and of (neuro)biology and biochemical sciences in general [START_REF] Kaila | pH and Brain Function[END_REF][START_REF] Zoulias | Hydrogen-Based Autonomous Power Systems[END_REF][START_REF] Dohle | [END_REF].

The common textbook example is centered on the water case because it clearly represents one of the most abundant substances on Earth and surely the most studied compound in condensed matter physics. In liquid water the process through which a H 2 O molecule dissociates, thus setting the acid or basic character of an aqueous solution, is known as protolysis and occurs according to the reaction

2H 2 O ⇋ OH -+ H 3 O + ,
in which, formally, a PT occurs between two water molecules.

The great importance played by this fundamental reaction is historically signed, among the others, by the fact that its basic peculiarities were elucidated before the discovery of the correct atomistic composition of the water molecule. In particular, in 1806 von Grotthuss [START_REF] Grotthuss | Sur la décomposition de l'eau et des corps qu'elle tient en dissolution à l'aide de l'électricité galvanique[END_REF], in an inspirational burst (trait de lumiére), develops a theory of water decomposition caused by galvanic electricity [START_REF] Cikierman | Et tu Grotthuss![END_REF]. Moreover, he suggested a primitive, but logically correct, diffusive model for the charges that he thought to compose the water molecule. The proposed original mechanism allows to account for the ultra-fast mobility of the proton H + which is about seven times than that of a sodium cation and five times than that of K + . Although the Grotthuss mechanism of charge transport has been proved also by many ab initio simulations [START_REF] Tuckerman | Ab initio Simulations of Water and Water Ions[END_REF][START_REF] Tuckerman | Ab initio Molecular Dynamics Simulation of the Solvation and Transport of H 3 O + and OH -Ions in Water[END_REF][START_REF] Tuckerman | Ab initio Molecular Dynamics Simulation of the Solvation and Transport of Hydronium and Hydroxyl Ions in Water[END_REF], the exact details of this fundamental process remain undisclosed and many models were presented in order to reproduce the available experimental evidences [START_REF] Agmon | The Grotthuss Mechanism[END_REF]. The tools provided by firstprinciples approaches are very suitable to tackle and treat chemical reactions which occur at this length-scale. Indeed by directly following the atomistic dynamical evolution it is possible to trace the relevant physical quantities and the local structural changes in assisting molecular dissociations and PT.

Atomic-scale computer experiments represent the state-of-the-art of computational science since they are capable to test very sophisticated theoriesi.e., quantum mechanics, quantum electrodynamics, etc. -via a direct comparison of their own predictions with real experimental results.

Since these kind of theories deal with many-body systems governed by difficult-/impossible-tosolve differential equations, computational physicists and mathematicians have developed smart and tricky algorithms and informatic tools in order to make concrete the possibility to simulate incredibly complex problems. However, the mere effort of computer sciences is not enough to account for the colorful multiplicity of ordinary condensed matter systems when treated in their most intimate details. Indeed, it has been necessary to develop many theoretical and mathematical apparata in order to make treatable the dynamical evolution of large-scale sets of molecules. Firstly, in order to make computationally solvable the electronic problem, Kohn and Sham (KS) [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF] gave one of the dominant contributions to the Density Functional Theory (DFT) by considering an auxiliary non-interacting (quantum) system sharing the same ground state density with the real interacting one. This way, an ensemble of fermions (electrons) can be (efficiently) computationally treated.

In addition to the just mentioned electron-electron part, the fundamental interactions that govern the nature of ordinary matter are ruled also by the proton-proton and the protonelectron terms. In most cases treating the behavior of a set of protons is order of magnitude less problematic than the electronic case. Indeed the quantum nature of the proton plays a central role only in a restricted number of phenomena and this allows us to use newtonian equations of motion for the nuclei; this classical nuclei approximation makes atoms computationally easily manageable. Moreover, Born and Oppenheimer [START_REF] Born | Quantum theory of the molecules[END_REF] have shown that there is no mixing between different electronic steady states caused by the interaction with the nuclei.

The KS scheme together with the last two introduced approximations for the nuclei define the Born-Oppenheimer Molecular Dynamics. This widely employed approach is the founder of the ab initio methods. Within this scheme the electronic problem has to be solved at each timestep and this leads to an increase of a factor greater than 1000 in the computational cost with respect to classical Molecular Dynamics. Historically this not negligible limitation of the ab initio simulations created more problems than it is doing nowadays because of the speed of growth of computing power. This circumstance led Roberto Car and Michele Parrinello in 1985 [START_REF] Car | Unified Approach for Molecular Dynamics and Density-Functional Theory[END_REF] to develop an incredibly smart ab initio method which overcomes the abovementioned time limitation. Their basic idea was to map the quantum-mechanical adiabatic time scale separation of fast electronic and slow nuclear motion into a classical-mechanical adiabatic energy-scale separation in the framework of dynamical system theory. This technique allows to tackle larger systems and for longer times than the previously employed first-principles methods, thus making feasible simulations of vastly different systems and phenomena. This approach has been thoroughly employed in most simulations carried out and here presented.

Coming back to our main subjects concerning molecular dissociation and PT, several experimental [START_REF] Meiboom | Nuclear Magnetic Resonance Study of the Proton Transfer in Water[END_REF][START_REF] Mohammed | Sequential Proton Transfer through Water Bridges in Acid-Base Reactions[END_REF][START_REF] Garczrek | Functional Waters in Intraprotein Proton Transfer Monitored by FTIR Difference Spectroscopy[END_REF][START_REF] Mafé | Electric Field-Assisted Proton Transfer and Water Dissociation at the Junction of a Fixed-Charge Bipolar Membrane[END_REF] and ab initio [6-8, 17, 18] studies were conducted in water samples. Despite the manifest importance of the liquid phase, investigating the PT in some of its crystalline structures can lead to new insights about this intricate phenomenon. Indeed the greater stability of the H-bonds and the increased atomic order of the solid phases with respect to the bulk liquid help to address more easily any detailed aspect in which molecular dissociation and PT can be decomposed. In this respect, ice phases are good candidates where to study both the static and the dynamical properties of the PT mechanism. At ambient conditions, the stable crystalline phase is ice I h , whose molecules are positionally arranged in an hexagonal lattice while being orientationally disordered (i.e., the proton sublattice is randomly distributed).

Starting from the standard ice crystalline structure it is possible to obtain its proton-ordered counterpart (ice XI) [START_REF] Tajima | Phase Transition in KOH-Doped Hexagonal Ice[END_REF][START_REF] Umemoto | Anomalous Pressure-Induced Transition(s) in Ice XI[END_REF], whose molecules are oriented along the c-axis (i.e., ferroelectric phase). The ordering of the proton sublattice obviously affects the properties of the material and, in particular, the electrical polarizability and the conductivity. In this respect ice can be described as a "protonic semiconductor". Indeed, ionic defects (i.e., hydronium and hydroxide) together with Bjerrum orientational defects are responsible for electrical transport processes in ice phases [START_REF] Petrenko | Physics of Ice[END_REF].

In the literature, some relevant efforts have been recently done to understand the PT mechanism, so as to fill the gap between theory and experiment [START_REF] Hassanali | Proton Transfer Through the Water Gossamer[END_REF][START_REF] Cwiklik | Hydroxide Trapped in the Interior of Ice: a Computational Study[END_REF][START_REF] Tuckerman | Structure and Dynamics of OH -(aq)[END_REF][START_REF] Dagrada | Quantum Monte Carlo Study of the Protonated Water Dimer[END_REF][START_REF] Bronstein | Quantum-driven Phase Transition in Ice Described via an Efficient Langevin Approach[END_REF]. One of the main difficulties is associated with the extremely rare nature of the protolysis phenomenon. In this respect, Saitta and coworkers [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF] have shown, in an ab initio study of water dissociation, that it is possible to stimulate the proton-transport process -and hence to investigate it in a systematic fashion -by applying an external electric field. In this thesis we extend this computational approach to crystalline phases of water. In fact, the peculiar symmetries of I h and XI ices allow a qualitative and quantitative investigation of the entropic and energetic contributions to the PT phenomenon.

As a logical consequence, in order to shed more light on this apparently elaborate phenomenon also by considering the ices mechanical response to the external field, we present in this thesis a series of further ab initio simulations performed to understand and exploit the symmetry differences between ice I h and ice XI. This aim is achieved by applying a static electric field along two different crystalline directions. We investigate several dynamical aspects of the H-bond structure during the proton migration process, which can be compared with the similar mechanism in liquid water. In this way it is possible to obtain a more exhaustive comprehension of the mechanical and electrical responses of ice I h and ice XI subject to the action of a static electrical perturbation.

Although the most important reactions which involve a PT between molecular entities occur certainly in aqueous environments, this phenomenon is so basic, fundamental, and vital for the occurrance of other related chemical reactions that its study should not be restricted to the water case. For this reason we extended our approach to the liquid methanol case. This way it is possible to highlight some of the relevant reaction coordinates which intervene in the complex interplay between intra-and inter-molecular interactions that assist the whole PT phenomenon in H-bonded systems. Single-donor H-bonds confer upon methanol an important role for the understanding of the PT mechanism and of the molecular dissociation process. Like water, liquid methanol displays a fast PT. In particular, it is believed that the PT process takes place along the H-bond structure via the Grotthuss mechanism [START_REF] Grotthuss | Sur la décomposition de l'eau et des corps qu'elle tient en dissolution à l'aide de l'électricité galvanique[END_REF], whose exact nature has been (and still is) a source of vivid controversies [START_REF] Agmon | Elementary Steps in Excited-State Proton Transfer[END_REF] and open debate.

Although its mere physico-chemical importance, this compound has a strong and multifaceted industrial involvement. As far as the proton conduction is concerned, methanol is used in apparatus for the production of electrical energy based on its fundamental molecular dissociation.

One example is provided by Direct Methanol Fuel Cell (DMFC) technology where electrical energy is produced by transferring protons [START_REF] Dohle | [END_REF]. Hence, a deeper understanding of the microscopic mechanism that leads to the molecular dissociation and/or to a net proton flow is of paramount importance to increase the efficiency of several ecological systems which produce energy from methanol.

In this thesis we present the first, to the best of our knowledge, ab initio Molecular Dynamics study of bulk liquid methanol under the effect of a static electric field and several analogies and differences with the corresponding liquid water case are highlighted.

As it is clear, the molecular dissociation event is extremely rare at ambient conditions. Indeed once an ionic pair is formed, the most likely circumstance is the ionic recombination. On average, the pH of ambient water implies that only about once, every 10 h, a randomly chosen water molecule will successfully dissociate spontaneously because of the local electric field fluctuations due to the molecular librations. A similar estimate can be obtained also for the methanol case since its dissociation constant is comparable with that of water. As already mentioned, in order to improve the statistical weight of the events in (almost) all the presented simulations, we applied a progressively increasing static and homogenoeus electric field, after several ps of Car-Parrinello Molecular Dynamics at ambient conditions. This way it was possible to stimulate the occurrance of molecular dissociations and to study the fundamental coupling between the molecules and the external field. Moreover, at slightly higher field strengths also the ionic conduction and the underlying PT mechanism could be analyzed.

Although our main interest was devoted to the detailed study of molecular dissociation and consequent PT in H-bonded systems, at moderate-to-high field strengths the chemical scenarios may be clearly more complex, as also recently proven in a prebiotic context [START_REF] Saitta | Experiments in Atomistic Computer Simulations[END_REF]. Of course, in the ice samples no chemical reactions were observed beyond the molecular dissociation/recombination but in the liquid methanol case the situation is deeply different. Indeed the atomistic reactivity of the carbon and oxygen atoms combined with intense field strengths leads to the onset of a multifaceted spectrum of "philicity" regions which are susceptible to nucleophilic or electrophilic attacks. Moreover, the electric field contributes to the cleveage of C-O, C-H, and O-H covalent bonds in the methanol molecules and subsequently to the formation of other combinations of these. Its entropic supply allows for the formation of species which are more or less complex than the methanol one such as formaldehyde (H 2 CO), methane (CH 4 ), or dimethyl ether (CH 3 OCH 3 ), and the evolution of these typical reactions can be predicted by exploiting many conceptual DFT tools such as the Fukui functions and indices [START_REF] Parr | Density Functional Approach to the Frontier-Electron Theory of Chemical Reactivity[END_REF][START_REF] Parr | Electronegativity: The Density Functional Viewpoint[END_REF].

In despite of the importance of the apparatus supplied by global and local (regio)selectivity parameters of modern conceptual DFT, one key aspect in analyzing any chemical reaction is represented by its energetics. The evaluation of the free energy suface (FES) can be tackled via disparate theoretical methods [START_REF] Voter | Accelerated Molecular Dynamics of Infrequent Events[END_REF][START_REF] Henkelman | A Climbing Image Nudged Elastic Band Method for Finding Saddle Points and Minimum Energy Paths[END_REF][START_REF] Torrie | Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling[END_REF] but the most reliable (and maybe used) ones are some metadynamics (MetD) approaches [START_REF] Laio | Escaping Free-Energy Minima[END_REF]. Basically, by introducing a history-dependent potential that fills the minima in the FES [START_REF] Laio | Escaping Free-Energy Minima[END_REF], it is possible to explore in principle every energetic landscape in the space spanned by few collective variables (CV) on which the free energy depends. In the last years different theoretical tools based on the MetD philosophy have been developed which circumvent the problem of choosing the correct CV (i.e., the irriducible set on which a specific FES relies upon). One of the most powerful is the path-CV one [START_REF] Branduardi; Gervasio | From A to B in Free Energy Space[END_REF] but it presents the drawback of knowing all the elementary steps via which a given reaction occurs before performing the MetD calculation. However, a new kind of path-CV MetD has been very recently presented [START_REF] Pietrucci | Formamide Reaction Network in Gas Phase and Solution via a Unified Theoretical Approach: Toward a Reconciliation of Different Prebiotic Scenarios[END_REF] that allows one to explore -and in some simple cases also to sample during the same simulation -the FES of in principle every kind of chemical reaction. The minimal a priori knowledge that this method requires are the simple approximated contact matrices of the atoms in the state of the reactants and in that of the products. From this perspective it represents an extremely powerful technique in exploring the FES of a chemical reaction. For more complex systemsi.e. FESs with many basins -the convergence of the MetD calculation is not so direct. Indeed, generally speaking, in principle the MetD exploration of the FES does not imply in any case its sampling! Hence once the "exploratory" FES has been obtained, other techniques must be employed in order to sample it. Umbrella sampling (US) [START_REF] Torrie | Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling[END_REF] is one of the most widely used methods for the determination of accurate free energies from atomistic simulations [START_REF] Torrie | Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling[END_REF][START_REF] Kumar | The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method[END_REF]. In this technique, independent simulations are carried out in different locations of the CV space with different biasing potentials in order to probe the underlying probability distributions. The just mentioned theoretical apparatus has been employed in order to accurately reproduce the energetics of some chemical reactions which were observed in previous Car-Parrinello Molecular Dynamics simulations in liquid methanol under the effect of a static electric field.

The present thesis is structured as follows. In Chapter 1 the main theoretical background thoroughly employed in this work will be presented. A less exhaustive appendix (Appendix A) has been also included in order to briefly accompany the reader in the conceptual DFT framework. In Chapter 2 we illustrate the results concerning the effects of applying moderate electric field intensities on three samples: ice I h , ice XI, and liquid methanol. In these systems instantaneous molecular dissociations occur and a net proton flow can be appreciated above determined field thresholds. At moderate/high field regimes several chemical reactions can be observed in the methanol sample and their study is carried out in Chapter 3 by invoking the theoretical apparatus shown in Appendix A and some marginal results presented in Appendix B.

Finally, some preliminary results on the energetics associated with one of the observed chemical reactions are shown at the end of Chapter 3. Beyond the Conclusions, the present thesis will close with the usual bibliographic references.

Chapter 1

Theoretical Background 1.1 Density Functional Theory

Hohenberg-Kohn Theorem

The fundamental theorem on which DFT is builded is the HK theorem [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF]. It represents the most basic of a number of existence theorems which ensure that stationary many-particle systems can be fully characterized by the ground-state density and intimately related quantities [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Engel | Density Functional Theory -An Advanced Course[END_REF]. This theorem, published in 1965, has been subsequently developed in order to include also systems at finite temperature [START_REF] Mermin | Thermal Properties of the Inhomogeneous Electron Gas[END_REF], systems under the action of magnetic field [START_REF] Von Barth | A Local Exchange-Correlation Potential for the Spin-Polarized Case[END_REF][START_REF] Pant | Theory of Inhomogeneous Magnetic Electron Gas[END_REF], and in general for all kinds of systems of (not only) physical domain.

Let us consider the standard Hamiltonian (in this paragraph in reduced units, and thoroughly in the thesis the convention 1 4πǫ 0 = 1 will be employed) of a stationary system of N interacting electrons under the action of the external potential v ext (r):

Ĥ = - 1 2 N i=1 ∇ i 2 + N i=1 v ext (r i ) + 1 2 N i=1 N j =i 1 |r i -r j | = T + Vext + Vee , (1.1)
where T is the kinetic operator, Vext is the external potential operator, and Vee is the operator which describes the electron-electron interaction. The first two terms on the right side of eq.

(1.1) are one-body operators whereas the third one is clearly a two-body operator. Let generalize the just introduced Hamiltonian by underspecializing the mutual electronic interaction by imposing Vext = Ŵ , where

Ŵ = 1 2 N i =j w(r i , r j ) = 1 2 σ,σ ′ =↑,↓ â † (rσ)â † (r ′ σ ′ )w(r, r ′ )â(r ′ σ ′ )â(rσ)d 3 rd 3 r ′ (1.2) whereas Vext = N i=1 v ext (r i ) = v ext (r)n(r)d 3 r . (1.3)
In eq. (1.2) r i is the position of the i-th particle, σ i is the z-projection of its spin, and â( †) (rσ) is the field operator which describes the creation/annihilation of one particle with spin σ at point r and which satisfies the general anticommutation rules. At this level is not necessary to specify the interaction w(r, r ′ ) in detail. The fundamental formalism of the whole DFT is independent of the form of w: it will be sufficient to assume that it is symmetric, w(r, r ′ ) = w(r ′ , r).

A : V → G .
Since the ground state density n 0 is

n 0 (r) = Ψ 0 |n 0 |Ψ 0 (1.4)
it can be defined a second relationship B between the set G and the set N of all the ground state densities resulting from some element of G via (1.4), i.e., N = {n

0 |n 0 (r) = Ψ 0 |n 0 |Ψ 0 , |Ψ 0 ∈ G }.
Ĥ|Ψ 0 = [ T + Vext + Ŵ ]|Ψ 0 = E 0 |Ψ 0 (1.5) Ĥ′ |Ψ 0 = [ T + V ′ ext + Ŵ ]|Ψ 0 = E ′ 0 |Ψ 0 . (1.6)
By subtracting side to side the two equations one obtains

[ Vext -V ′ ext ]|Ψ 0 = [E 0 -E ′ 0 ]|Ψ 0 . (1.7)
If we rewrite the eq. (1.7) in the first quantized formalism rather than in that of the second, it is possible to divide by the ground state wavefunction thanks to the multiplicative character (i.e., local) of the potential and thus one can obtain:

N i=1 [v ext (r i ) -v ′ ext (r i )] = E 0 -E ′ 0 . (1.8)
For all the r i such that the wavefunction does not vanish, let keep fixed N -1 of them in eq.

(1.8), whereas let vary the remainder position. This process leads to an evident contradiction since at the right-hand side, by construction of the set V , we have a quantity which cannot be a constant, whereas at the left-hand side we have actually a constant. As a consequence the map A is injective. Now let us turn to the map B. Its definition implies that it is impossible to relate |Ψ 0 to two distinct n 0 and n ′ 0 . This way, in order to prove the injectivity of B, one has to show that two different |Ψ 0 ∈ G cannot result in the same n 0 . The proof relies again on a reductio ad absurdum: let us assume that |Ψ 0 and |Ψ ′ 0 ∈ G are linked to the same n 0 . From the Ritz variational principle for the ground state energy one has

E 0 = Ψ 0 | Ĥ|Ψ 0 < Ψ ′ 0 | Ĥ|Ψ ′ 0 , (1.9)
where Ĥ is the unique Hamiltonian leading to |Ψ 0 (due to the uniqueness of A) and the strict inequality is due to the non-degeneracy of |Ψ 0 and |Ψ ′ 0 . The right-hand side of (1.9) can be evaluated by adding and subtracting the unique potential V ′ ext corresponding to |Ψ ′ 0 :

E 0 < E ′ 0 + Ψ ′ 0 | Vext -V ′ ext |Ψ ′ 0 . (1.10)
Because of the local character of the potential and due to the (1.3) one obtains

E 0 < E ′ 0 + n 0 (r)[v ext (r) -v ′ ext (r)]d 3 r . (1.11)
However, the same argumentation could be further conducted by exchanging the apex which distinguishes the distinct potentials leading in such a case to

E ′ 0 < E 0 + n 0 (r)[v ′ ext (r) -v ext (r)]d 3 r . (1.12)
Upon addition of (1.11) and (1.12) one arrives to a clear contradiction:

E 0 + E ′ 0 < E ′ 0 + E 0 . (1.13)
Is is possible to conclude that the map B is unique: it exists a one-to-one correspondence between |Ψ 0 and n 0 . Q.E .D.

Physical Insights

The HK theorem has a series of physical implications of primary importance. Indeed, due to the uniqueness of the relationships between v ext , |Ψ 0 , and n 0 , the ground state is a unique functional of the ground state density and it can be labeled as |Ψ 0 = |Ψ 0 [n 0 ] by inverting the map B (i.e., B -1 ). It has been noticed that no assumptions were made on the specific form of

v ext in constructing |Ψ 0 [n 0 ] from n 0 ; |Ψ 0 [n 0 ]
has the same functional form for every kind of many-particle system with the same interaction w. |Ψ 0 [n 0 ] is then a universal functional from which it is possible to determine every physical observable,

O[n] = Ψ[n]| Ô|Ψ[n] , (1.14)
such as the ground state energy

E[n] = Ψ[n]| Ĥ|Ψ[n] = F [n] + v ext (r)n(r)d 3 r , ( 1.15)
where

F [n] = Ψ[n]| T + Ŵ |Ψ[n] (1.16)
is the universal part. It exists a minimum principle for E[n], according to which

E[n 0 ] < E[n ′ 0 ] ⇐⇒ E 0 = min n∈G E[n] .
(1.17)

Hence the E[n] domain is restricted to the densities belonging to N ; only the ground state densities stemming from the Schrödinger equation are legitimate and not all the arbitrarily functions of r which are non-negative and with finite norm.

From a strictly physical perspective of course the particle react to every arbitrarily small change of v ext , such that their total energy can be minimized, and this response is unique. Because of the local character of the potential, no changes will occur in v ext without a modification of the electronic wave function and thus of the density distribution.

The Levy approach and the Kohn-Sham scheme

Up to now we do not faced the problem of v-representability or of the N -representability because the HK theorem, in despite of its conceptual and historical importance, is not the best starting point for this topic. A very useful approach is that of the constrained-search introduced by Levy [START_REF] Levy | Universal Variational Functionals of Electron Densities, First Order Density Matrices, and Natural Spin-Orbitals and Solution of the v-Representability Problem[END_REF][START_REF] Levy | Electron Densities in Search of Hamiltonians[END_REF] and further developed by Lieb [45]. Eq. (1.17) can be explicitly derived from the following variational problem

δ    Ψ| Ĥ|Ψ Ψ|Ψ    = 0 . (1.18)
This equation highlights that the ground state energy can be found by minimizing Ψ| Ĥ|Ψ over all the N-particle antisymmetric normalized wavefunctions:

E = min Ψ Ψ| Ĥ|Ψ . (1.19)
Let us divide the minimization process of this expression in two steps.

1. Let us consider all the wavefunctions which produce a given density n(r) and let us minimize the Hamiltonian with respect to these eigenfunctions:

min Ψ→n Ψ| Ĥ|Ψ = min Ψ→n Ψ| T + Vee |Ψ + v ext (r)n(r)d 3 r (1.20)
where we replaced the electron-electron interaction Vee and where we exploited the fact that all the |Ψ which leads to the same n(r), leads also to the same Ψ| Vext |Ψ . Let us now redefine the universal functional of (1.16) in the restricted domain of all the eigenfunctions |Ψ min n that minimize the value of T + Vee for a given n(r).

2. Finally, let us minimize over all the N-electrons densities n(r):

E = min n E[n] = min n F [n] + v ext (r)n(r)d 3 r (1.21)
where v ext (r) is kept fixed. The density that minimizes the term in the parenthesis is the ground state density.

The constraint of keeping fixed the number of electrons N can be managed by means of the Lagrange multiplier µ:

δ F [n] + v ext (r)n(r)d 3 r -µ n(r)d 3 r = 0 (1.22)
which is equivalent to the Euler equation:

δF δn(r) + v ext (r) = µ , (1.23)
in which µ will be adjusted until the equation for the (electronic) mass conservation (i.e.,

σ n σ (r)d 3 r = N ) is satisfied. The Euler equation (1.23) shows that the external potential v ext (r) is uniquely determined by the ground state density (or by one of them in case of degeneracy). The functional F [n] of the expression (1.16) is defined for all the densities n(r) which are "N -representable"; in other words, for the densities stemming from a N-electrons antisymmetric wavefunction. The functional derivative δF [n] δn(r) is defined by the Euler equation for all the densities n(r) which are "v-representable"; in other words, for the densitites stemming from ground state N-electrons antisymmetric wavefunctions corresponding to a given external potential v ext (r).

The Levy approach is also the most natural method which allows one to introduce the Kohn-Sham (KS) scheme [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]. The basic idea arises from the necessity of solving the static and the dynamics of a set of interacting electrons. Kohn and Sham considered an auxiliary electronic system (i.e., the KS system) which was characterized by the same electron density of the "real" one but without interactions between the fermions. This mathematical artifact allows one to remove the operator associated with the mutual electronic interactioni.e., Vee = 0 -by displacing the focus on the fictitious potential v s (r) (KS potential) that produces the same density and that satisfies the Euler equation (1.23) in this fashion:

δT s δn(r) + v s (r) = µ (1.24)
where

T s [n] = min Ψ→n Ψ| T |Ψ = Φ min n | T |Φ min n . (1.25)
The minimizing wavefunction |Φ min n for a given density should be a non-interacting one in the sense that it should be obtained as a product of one-particle states such as a single Slater determinant or a linear combination of them, correspondingly to an external potential operator

Vs . In eq. (1.24), the KS potential v s (r) is a functional of n(r) and, if it was present a difference between µ and µ s (i.e., the chemical potentials of the interacting and non-interacting systems, respectively), it would be adsorbed in v s (r).

Let us define now the exchange-correlation (XC) energy E xc [n]:

F [n] = T s [n] + U [n] + E xc [n] (1.
v s (r) = v ext (r) + δU [n] δn(r) + δE xc [n] δn(r) . (1.27)
It has to be noticed that by means of the functional derivative process the KS potential can be set in the form

v s (r) = v ext (r) + v H [n 0 ](r) + v xc [n 0 ](r) , (1.28)
where the external potential, the Hartree contribution (i.e., purely electrostatic repulsion), and the XC potential are present. Clearly, the unknown potential v s is by itself a density functional, thus the Schrödinger equations result to be non linear; the solutions |ψ i determine the density n 0 by means of the relationship

n 0σ (r) = σ=↑,↓ i Θ(µ -E iσ )|ψ i (rσ)| 2 , (1.29)
where i is the set of quantum numbers (regardless of the spin component) and Θ is the Heaviside step function which ensures that all the KS spin orbitals with E iσ < µ are singly occupied because of the Pauli esclusion principle and that those with E iσ > µ are empty. Hence, the previous equation determines the KS potential v s via eq. (1.28). The KS equations result from the substitution of eq. (1.28) in the one-particle Schrödinger equation,

- 2 2m ∇ 2 + v ext (r) + v H [n](r) + v xc [n](r) ψ i (rσ) = ε i ψ i (rσ) (1.30)
and they has to be solved self-consistently. The starting point is a suitable choise of a given test density n (1) which allows to build the initial KS potential (1) ](r) + v xc [n (1) ](r) (1.31) provided that the functional derivative of E xc [n] is known. This KS potential will be employed in order to solve eq. (1.30) which produces some wavefunctions ψ

v (1) s (r) = v ext (r) + v H [n
(2) i which in turn lead to an "improved" density by means of eq. (1.29). Density n (2) may in turn be used to obtain other potentials, v H [n (2) ](r) e v xc [n (2) ](r), which better fit the real system. This iterative procedure has to be continued until the difference between two densities obtained in two successive iterations is lower than a previously chosen threshold of accuracy.

The exchange and correlation functional

The XC functional E xc [n] is the sum of two distinct terms:

E xc [n] = E x [n] + E c [n] , (1.32)
where

E x [n] = Ψ min n | Vee |Ψ min n -U [n] . (1.33)
If Ψ min n is a single Slater determinant, eq. (1.33) is the usual Hartee-Fock applied to the KS orbitals. As an example, it differs from the HF exchange energy at the same way in which the KS orbitals or T s [n] differ from the HF orbitals and from the HF kinetic energy, respectively.

The correlation energy can be defined as .34) In this latter Ψ min n is the eigenfunction that produces the density n and that minimizes T + Vee , whereas Φ min n is the eigenfunction which corresponds to the same density but that minimizes T . It is clear that, by the definition of the universal functional F [n], the correlation energy is non-positive.

E c [n] = F [n] -T s [n] + U [n] + E x [n] = Ψ min n | T + Vee |Ψ min n -Φ min n | T + Vee |Φ min n . ( 1
The XC functional is then necessary in order to include all that is not involved, within the KS scheme, in the kinetic term T s , in the mutual electronic interaction, and in the external potential, in order to determine the ground state energy. Basically if the exact dependence of the XC functional on the densities n ↑ and n ↓ was known, the KS equations would predict the exact ground state energy and the exact spin densities of a many-electron system. The forces acting on the nuclei and their equilibrium positions could be then evaluated by means of the spatial gradient of this energy with the sign changed.

In practice E xc has to be approximated. For a long while the most used approximation was the Local Spin Density Approximation (LSD) [START_REF] Von Barth | A Local Exchange-Correlation Potential for the Spin-Polarized Case[END_REF] or, in the spin-unpolarized case the Local Density Approximation (LDA) [START_REF] Kohn | Self-consistent equations including exchange and correlation effects[END_REF]; by considering the spin densities, one has

E LSD xc [n ↑ , n ↓ ] = n(r)e xc (n ↑ (r), n ↓ (r))d 3 r , (1.35)
where e xc (n ↑ (r), n ↓ (r)) is the XC energy per particle referred to an electronic gas with uniform spin densities n ↑ and n ↓ and which is known. By treating the electronic density as homogeneous, this approximation cannot take into account the high complexity of the energy functional dependence on the densities. More recently, approximations such as the Generalized Gradient Approximations (GGA) have been developed which consider also the functional dependence on the density gradient:

E GGA xc [n ↑ , n ↓ ] = f (n ↑ , n ↓ , ∇n ↑ , ∇n ↓ )d 3 r . (1.36)
Even if E xc is often a relatively small fraction of the total energy of an atom or a molecule, its contribution is comparable to the bond energy or to the atomization one. Without this term the atoms would be in the most favored cases only weakly bonded. It is clear that a central issue of the whole DFT framework is to have more and more accurate approximations, such as those hybrid [START_REF] Becke | Density-functional Thermochemistry. 3. The Role of Exact Exchange[END_REF][START_REF] Becke | A New Mixing of Hartree-Fock and Local Density Functional Theories[END_REF] that include to some extent Hartree-Fock (HF) exchange or full exact exchange (EXX) [START_REF] Grimm | Restricted Open-shell Kohn-Sham Theory for π-π * Transitions. I. Polyenes, Cyanines and Protonated Imines[END_REF] in addition to the standard GGA expression. Another schemes which go beyond GGA are the meta-GGA (MGGA) [START_REF] Tao | Climbing the Density Functional Ladder: Nonempirical Meta-Generalized Gradient Approximation Designed for Molecules and Solids[END_REF] (which include higher-order powers of the gradient) and the hyper-GGA (HGGA) [START_REF] Perdew | Prescription for the Design and Selection of Density Functional Approximations: More Constraint Satisfaction with Fewer Fits[END_REF] (which uses, in addition to the characteristics of the MGGAs, the exact exchange energy density).

In this thesis we will apply mainly the Perdew-Burke-Ernzerhof (PBE) [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] and secondly the Becke-Lee-Yang-Parr (BLYP) [START_REF] Becke | Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior[END_REF] functionals, both belonging to the GGA class.

PBE

Let us at first remember the exchange part of the energy density of the homogeneous electron gas (HEG):

ǫ X HEG [n] = - 3 4 3 π 1/3 n 1/3 , (1.37)
whereas for the correlation part one can split it into two terms depending on the density regime (i.e., its solutions are known only in the low-and high-density limits).

The analytical expression of the correlation part of the so-called PBE functional reads as

E C P BE [n] = n(ǫ C HEG (n) + H C [n, t])dr , t = ∇n 2k s n (1.38)
where t is a dimensionless density gradient, 1/k s is the Thomas-Fermi screening wave-length, and

H C (n, t) = a 2 m 2 ln 1 + a 1 a 2 t 2 1 + At 2 1 + At 2 + A 2 t 4 , ( 1.39)
where

A = a 1 a 2 1 exp (-ǫ C HEG /(a 2 m/ 2 )) -1 (1.40)
with a 1 = 0.066725 and a 2 = (1 -ln 2)/π 2 .

For the exchange part one has:

E X P BE [n] = n(r)ǫ X HEG [n]F X (s)dr , s = 2(3π 2 /2) 1/3 √ r s t , (1.41)
where s is a dimensionless gradient and

F X (s) = 1 + b 1 - b 1 1 + b 2 s 2 b 1 (1.42)
with b 1 = 0.804 and b 2 = a 2 π 2 /3.

BLYP

Becke proposed to correct the exchange density from LDA (eq. 1.37) in the following manner

E X B [n] = nǫ X HEG [n] 1 + 4 3 π 3 1/3 a 1 x 2 1 + 6a 1 x sinh -1 x dr , x = |∇n| n 4/3 , (1.43)
where a 1 = 0.0042 is obtained via a least square fit on the exact HF exchange of noble gases.

The well-known LYP correlation functional has been obtained from a semplification and a reformulation of the Colle-Salvetti formula in terms of n and ∇n only:

E C LY P [n] = - b 1 1 + b 4 n -1/3 n + b 2 n -2/3 C F n 5/3 -2t w + t w 9 + ∇ 2 n 18 e -b 3 n -1/3 dr , (1.44)
where

C F = 3 10 (3π) 2/3 , t w = 1 8 |∇n| 2 n -∇ 2 n .

Plane waves and pseudopotentials

The total energy of a solid can be evaluated through a direct calculation of the KS wavefunctions. This is very useful since the periodicity of the problem can be clearly exploited. For finite systems, such as molecules, clusters and so on, it is always possible to take advantage of the plane waves mathematical properties in a super-cell approach. By considering the infinite replication of the super-cell it is quite straightforward to solve the KS equations in the Fourier space. Moreover, in many cases one is interested to small variations of the total energy associated with the valence electrons rearrangement regardless of the core contribution. Indeed inner electrons are often chemically inert and do not significantly contribute in bonding. This led to the introduction of pseudopotentials felt by the valence electrons and which are generated by the joint effect of the nuclear charge and the inner electrons.

Within a super-cell approach and by making use of a plane waves basis set, it is possible to divide the solutions of the Schrödinger equation Ĥ|Ψ n = E n |Ψ n in two sets, |Ψ c and |Ψ v , i.e., for the inner electrons and the valence ones, respectively (n = c, v). The valence orbitals can be written as the sum of a smooth function |φ v , the pseudofunction, with an oscillating function which results from the orthogonalization of the valence to the inner orbitals

|Ψ v = |φ v + c α cv |Ψ c , (1.46)
where α cv = -Ψ c |φ v . The Schrödinger equation for the smooth orbital |φ v leads to

Ĥ|φ v = E v |φ v + c (E c -E v )|Ψ c Ψ c |φ v . (1.47)
This latter highlights the fact that the states |φ v satisfy a Schrödinger-like equation with an energy-dependent pseudo-Hamiltonian

ĤP (E) = Ĥ - c (E c -E v )|Ψ c Ψ c | . (1.48)
As a consequence

ωP (E) = v - c (E c -E)|Ψ c Ψ c | , (1.49)
where v is the true potential in which valence electrons move. Far away the core region ωP becomes v. The valence electrons are subject to an effective potential which is the combination of the screening of the nuclear charge due to the core electrons, the Pauli repulsion, and the XC effects between the valence and core electrons. Generally the pseudopotential can be cast into the form

ω(r, r ′ ) = l l m=-l Y lm * (r)ω l (r, r ′ )Y lm (r ′ ) , (1.50)
where Y lm are the spherical harmonics.

In a plane wave basis set calculation and with pseudopotentials, the external potential v ext is simply the sum of all pseudopotentials of all atoms in the system. If atom α is located in the unit (super-)cell at r α and its pseudopotential is

ω α (r, r ′ ), v ext becomes ω(r, r ′ ) = j,α ω α (r -R j -r α , r ′ -R j -r α ) , (1.51)
where R j are the lattice vectors.

The KS wavefunctions ψ i (r), according to the Bloch's theorem will be

ψ k,n (r) = e ik•r G c k,n (G)e iG•r , ( 1.52)
where k is the wavevector, n the band index, and G are the reciprocal lattice vectors. From this latter expression one can obtain the density and its Fourier transform

n(G) = k,n G ′ f (ǫ k,n )c k,n * (G ′ -G)c k,n (G ′ ) . (1.53)
The kinetic energy is thus rewritten as

T = 1 2 k,n G f (ǫ k,n )|c k,n (G)| 2 |k + G| 2 , ( 1.54)
and the Hartree energy

E H = Ω 2 G v H (G)n(G) , (1.55)
where Ω is the unit cell volume. It has to be noticed that one way to obtain the Hartree potential is through the respective Poisson equation.

The electron-nuclei interaction energy can be then expressed as

E e-n = k,n G,G ′ f (ǫ k,n )c * k,n (G)c k,n (G ′ )ω(k + G, k + G ′ ) . (1.56)
E e-n , the Hartree potential, and the nucleus-nucleus interaction energy diverge at G = 0 but their sum is a constant value equal to E rep + E Ewald , if the system is electrically neutral [START_REF] Ihm | Momentum Space Formalism for the Total Energy of Solids[END_REF][START_REF] Ihm | Total Energy Calculations in Solid State Physics[END_REF][START_REF] Pickett | Pseudopotential Methods in Condensed Matter Applications[END_REF],

where

E rep = Z tot (1/Ω) α Λ α and Λ α = 1 Ω v α,local (r) + Z α r . (1.57)
Z α is the electric charge of nuclei α and v α,local is the local part of the pseudopotential of atom α. By separating the short-and the long-range non-divergent parts of the nuclei-nuclei interaction, E Ewald , one obtains

E Ewald = 1 2 α,α ′ Z α Γ α,α ′ Z α ′ , ( 1.58)
where Γ α,α ′ is the typical term appearing in the Ewald summations over the reciprocal lattice vectors and on the specific atomic positions and which does not depend on the density.

Finally, the total energy in the Fourier space is given by

E tot = T + E ′ Hartree + E ′ e-n + E xc + E Ewald + E rep , (1.59)
with the terms G,G ′ = 0 excluded from the Hartree and pseudopotential contributions. Hence the KS equations become

G ′ ĤG,G ′ (k)c k,n (G ′ ) = ǫ k,n c k,n (G) , ( 1.60)
where

ĤG,G ′ (k) = 1 2 |k + G| 2 δ G,G ′ + ω(k + G, k + G ′ ) + v H (G -G ′ ) + v xc (G -G ′ ) , (1.61)
and are solved by diagonalizing the Hamiltonian.

The adoption of plane waves basis sets even for systems that are not ordered can be justified in view of the important consequences that its use has in ab initio molecular dynamics (AIMD) simulations. Indeed as shown before, being plane waves a complete and orthonormal set of functions, they can used in order to expand the orbitals including only those G vectors that satisfy the periodic boundary conditions (PBC). Moreover, by introducing this originless basis set one avoid the onset of spurious forcesi.e., the Pulay forces -which arises from atom-centered (or localized) basis sets such as, e.g., Gaussians basis sets. This advantage, which occurs even within a finite basis and a restricted number of plane waves, leads in AIMD simulations to an extremely lower computational demand in particular in evaluating forces and stress tensors.

The simple and only way to improve the quality of the basis set is to increase the so-called "energy cutoff". Indeed the KS potential converges very rapidly with increasing of the modulus of G. Hence, at each k-point, only G vectors with a kinetic energy lower than a given maximum cutoff

1 2 |k + G| 2 ≤ E cut (1.62)
will be included in the basis.

Finally, plane waves seems really to be one of the most natural way of implementing basis sets in numerical calculations. Indeed, being defined very easily in the Fourier space, all the related calculations can be made here and the derivatives of real space becomes multiplications. This leads again to a computational saving and the results in the wavevector space can be efficiently mapped into the real one via Fast Fourier Transforms (FFTs).

Ab initio simulations

Up to now, we considered the nuclei as static entities since they are involved in pure DFT calculations only as charged particles in the external potential operator. In order to obtain the dynamical evolution of a set of atoms or molecules one has to derive specific equations of motion that take into account all the relevant degrees of freedom.

The Hamiltonian of any real system cannot transcend on the nuclear variables; in its more general structure it can be written as

Ĥ = - P I=1 2 2M I ∇ 2 I - N i=1 2 2m ∇ 2 i + e 2 2 P I =J Z I Z J |R I -R J | + e 2 2 N i =j 1 |r i -r j | -e 2 P I=1 N i=1 Z I |R I -r i | , ( 1.63)
where {R I } is the set of the P nuclear coordinates and {r i } is that of the N electronic ones, while Z I and M I are the nuclear charges and masses, respectively. An analytical solution to the Schrödinger equation corresponding to the Hamiltonian (2.63) can be achieved only in very few cases. Indeed, because of the correlations between all the involved particles one has to deal with 3(N + P ) coupled degrees of freedom. In order to simplify the electronic part of the problem one can introduce the just treated pseudopotentials and, for what concern the nuclear part, two approximations were proposed in first-principles calculations: the adiabatic and the classical nuclei approximations.

The first one stems from a perturbative analysis of the stationary Hamiltonian (2.63) with respect to the mass ratio κ = (m/M ) 1/4 and by considering a new nuclear coordinates set defined as R = R 0 + κu, where u represents the displacement of the nuclei from the equilibrium position R 0 . Under suitable conditions, Born and Oppenheimer proved the practical abscence of mixing between distinct electronic states caused by the interaction with the nuclei [START_REF] Born | Quantum theory of the molecules[END_REF]. In other words, while the nuclei follow their own dynamics, the electrons istantaneously adjust their total wavefunction according to the total nuclear wavefunction: these states are adiabatic.

Thus, under proper external conditions it is legitimate to separate the electronic and the nuclear variables. It has been shown that the correction to the evaluation of the eigenvalues is proportional to κ 4 = m/M , which in the worst cases is lower than the 0.5% [START_REF] Landy | The Adiabatic Approximation[END_REF]. Neglecting these corrective terms constitutes the Born-Oppenheimer approximation.

The classical nuclei approximation is based on the observation that the de Broglie wavelength of most of the nuclei is comparable with the interatomic distances only at temperatures of the order of few units of Kelvin. Hence the quantum coherence effects associated with the nuclei are extremely relevant only in these regimes (or almost only in these cases!). This circumstance allows one to treat the nuclei as classical objects and to compute the relative forces via some variants of the Hellmann-Feynman theorem [START_REF] Hellmann | Einfürung in die Quantenchimie[END_REF][START_REF] Feynman | Forces in molecules[END_REF]. The set of newtonian equations of motion for the nuclear degrees of freedom is

M I d 2 R I (t) dt 2 = - ∂ ∂R I Ψ(R)| ĥe |Ψ(R) - ∂V nn (R) ∂R I , ( 1.64)
where ĥe = -

2 2m N i=1 ∇ 2 i + e 2 2 N i =j 1 |r i -r j | + P I=1 N i=1 v ext (r i -R I ) (1.65)
and

V nn (R) = e 2 2 P I =J Z I Z J |R I -R J | , ( 1.66)
can be numerically integrated in order to obtain physical trajectories in the phase space. Within this scheme the electronic problem has to be solved at each integration step of the molecular dynamics (MD), tipically leading to an increase of the computational cost of a factor of 1000 or more with respect to the classical simulations. This kind of simulation is mandatory when the chemistry of the system plays a relevant role, such as in chemical reactions, during which the breaking and the reforming processes of covalent bonds take place.

A self-consistent first-principles MD approach can be obtained by solving the electronic problem within the KS scheme:

Ψ(R)| ĥe (R)|Ψ(R) = E KS [n](R) = T s [n] + 1 2 n(r)n(r ′ ) |r -r ′ | d 3 rd 3 r ′ + E XC [n] + P I=1 n(r)v ext (r -R I )d 3 r + 1 2 P I =J Z I Z J |R I -R J | , ( 1.67)
and the force on the nuclear coordinates can be obtained by simple derivation:

F I = - ∂E KS [n](R) ∂R I = -n(r) ∂v ext (r -R I ) ∂R I d 3 r + P I =J Z I Z J R I -R J |R I -R J | 3 . (1.68)
The standard procedure before the Car and Parrinello seminal article [START_REF] Car | Unified Approach for Molecular Dynamics and Density-Functional Theory[END_REF], was that of minimizing the electronic structure at each step of the MD simulation. This method is known as Born-Oppenheimer MD (BOMD) [START_REF] Bendt | Simultaneous Relaxation of Nuclear Geometries and Electronic Charge Densities in Electronic Structure Theories[END_REF] and several calculations presented in this thesis have been performed also within this approach. But at the time in which this method has been set up, in 1983, requiring the self-consistency of the electronic wavefunctions at each step represented the limiting step of the size and/or of the time length of the simulated sample.

Car-Parrinello Molecular Dynamics

The massive computational implementation of DFT started historically since 1985 thanks to the Car-Parrinello (CP) method [START_REF] Car | Unified Approach for Molecular Dynamics and Density-Functional Theory[END_REF]. Indeed it was presented an alternative approach to the BOMD by considering the possibility to maintain the electronic density sufficiently close to the adiabatic density during distinct MD steps. By introducing a fictitious conservative term at the second order of the dynamics of the KS orbitals, this purpose has been reached. The CP molecular dynamics (CPMD) scheme can be described by a dynamical system represented by the following Lagrangian:

L CP = 1 2 P I=1 M I Ṙ2 I + µ N i=1 f i | ψi (r)| 2 d 3 r -E KS [ψ i (r)](R) + N i=1 f i N j=1 Λ ij ψ * i (r)ψ i (r)d 3 r -δ ij , (1.69)
where the KS energy is given by

E KS [ψ i (r)](R) = N i=1 f i ψ * i (r) ĥ1 (R, r)ψ i (r)d 3 r + 1 2 n(r)n(r ′ ) |r -r ′ | d 3 rd 3 r ′ + E XC [n] + 1 2 P I =J Z I Z J |R I -R J | . (1.70)
The operator ĥ1 is the one-electron Hamiltonian in which are present the free kinetic term and the external potential, whereas the electronic density in terms of the dynamical KS orbitals and of their occupation numbers can be written as

n(r) = N i=1 f i ψ * i (r)ψ i (r)d 3 r . (1.71)
The first term in the CP Lagrangian (2.69) represents the nuclear kinetic energy while the last one guarantees the orthonormality, which is necessary to describe electronic densities arising from a Slater determinant of the KS orbitals at each simulation time step. The Lagrange multipliers Λ ij should be determined in such a way that this condition is satisfied. From a computational point of view, there exist many orthonormalization algorithms. Some of them do not guarantee the energy conservationi.e., the Gram-Schmidt algorithm -and are exploitable only in the earliest minimization steps of the energy. On the other hand, within a CP dynamics it exists only one orthonormalization procedure which is consistent with the equations of motion [START_REF] Pastore | Theory of Ab initio Molecular Dynamics Calculations[END_REF].

The second term in the CP Lagrangian (2.69) represents (with the orthonormalization constraints) the most important innovation introduced. It corresponds indeed to a fictitious kinetic energy

K e [ψ i (r)] = µ N i=1 f i | ψi (r)| 2 d 3 r (1.72)
associated with the dynamical evolution of the KS orbitals. K e is simply a convenient term that allows the KS orbitals evolution independently of (but coupled with) the nuclear degrees of freedom. K e is a measure of the deviations from the exact BO surface and in the limit µ → 0 one recovers the BOMD scheme. The fictitious mass, which is also known as adiabaticity or inertia parameter, assigned to the orbital degrees of freedom has the units of an energy times a squared time only for dimensionality reasons.

The equations of motion stemming from (2.69) are

d dt ∂L CP ∂ ṘI = - ∂L CP ∂R I , (1.73) d dt δL CP δ ψ * i (r) = - δL CP δψ * i (r)
.

(1.74)

The second equation involves functional derivatives because the wavefunctions are continuous scalar fields. However, in practice the wavefunctions are expressed with respect to a discrete basis or in terms of plane waves components leading to partial derivatives. This way one obtain the CP equations of motion

M I RI = - ∂E KS [ψ i (r)](R) ∂R I , (1.75) µ ψi (r, t) = - δE KS [ψ i (r)](R) δψ * i (r) + N j=1 Λ ij ψ j (r, t) = -ĤKS ψ i (r, t) + N j=1 Λ ij ψ j (r, t) . (1.76)
The introduction of the electronic inertia parameter µ leads obviously to a twofold temperature concept in CP simulations. Indeed the instantaneous physical temperature, which is proportional to I M I ṘI , associated with the nuclear motion has now the fictitious electronic counterpart ∝ i µ ψ i |ψ i . A central point is related to the decoupling of the two energy scales because during a CP simulation no energy transfer must be observed between the nuclear and the electronic subsystems. By maintaining the electrons sufficiently "cold" and the nuclei "hot" enough (but not too much!) it is possible to prevent a non-adiabatic time evolution and this aim is in practice achieved by accurately checking the relative power spectra. As an example, by considering the Fourier transform of the statistically averaged velocity autocorrelation function of the classical fields representing the orbitals

f (ω) = ∞ 0 cos (ωt) i ψi ; t| ψi ; 0 dt (1.77)
and by comparing it to the highest frequancy phonon mode ω max n , one can establish if the adiabaticity will be fulfilled. Indeed if no overlap is present between the two spectra, no energy transfer can occur between the two subsystems according to nonlinear classical dynamics theory.

In practice a harmonic analysis should be sufficient in evaluating this particular issue [START_REF] Pastore | Theory of Ab initio Molecular Dynamics Calculations[END_REF]. It has been shown [START_REF] Pastore | Theory of Ab initio Molecular Dynamics Calculations[END_REF] that the frequency spectrum of the classical orbitals close to the ground state is

ω ij =   2(ǫ i -ǫ j ) µ   1/2 , ( 1.78)
where ǫ i and ǫ j are the eigenvalues of unoccupied and occupied orbitals, respectively. By following this simple harmonic approximation [START_REF] Pastore | Theory of Ab initio Molecular Dynamics Calculations[END_REF], one can determine the lowest electronic "vibrational" frequency:

ω min e ∝   E gap µ   1/2 . (1.79)
This latter shows that by taking greater and greater values of µ, the minimum electronic frequency shifts to lower values increasing the overlap probability with the phonon spectrum.

Hence, in order to guarantee the adiabaticity one has to choose a value of µ small enough to have the biggest difference ω min e -ω max n as possible. However, the fictitious mass parameter has other responsabilities. For instance, decreasing µ not only shifts the lower frequency at higher values but also stretches the whole frequency spectrum because of the following relation

ω max e ∝   E cut µ   1/2 , (1.80)
where E cut is the largest kinetic energy in a plane waves expansion. Moreover a great limitation in arbitrarily decreasing the value of µ is due to its relationship with the allowed maximum time step

∆t max ∝   µ E cut   1/2 . (1.81)
Hence, depending on the specific system under investigation and by playing with µ and E cut one can set up the desired CP simulation under adiabatic conditions.

Of course, the insertion of µ has also several drawbacks. The primary effects is clearly due to the fact that the electronic subsystem has a mass! This means that the electronic degrees of freedom exert some inertia on the nuclei making them effectively heavier and this circumstance manifests itself sometimes mainly in thermodynamic and structural properties. Fortunately, by means of a perturbative treatment of the CP single-particle wavefunctions with respect to the corresponding BO ones [START_REF] Marx | Ab Initio Molecular Dynamics[END_REF][START_REF] Tangney | How well do Car-Parrinello Simulations Reproduce the Born-Oppenheimer Surface? Theory and Examples[END_REF], one recovers not only suitable relationships between the respective forces, but also renormalization values for the nuclear masses and for the relative phonon frequency scaling; this a way, a proper temperature scaling can be adapted.

Static electric fields in ab initio simulations

The theoretical development and the implementation of electric fields in ab initio codes is all other than trivial. A wide literature exists in this field [START_REF] Nunes | Real-space Approach to Calculation of Electric Polarization and Dielectric Constants[END_REF][START_REF] Nunes | Berry-phase Treatment of the Homogeneous Electric Field Perturbation in Insulators[END_REF][START_REF] Umari | Ab initio Molecular Dynamics in a Finite Homogeneous Electric Field[END_REF][START_REF] Resta | Quantum-mechanical Position Operator in Extended Systems[END_REF][START_REF] Resta | Chapter of Physics of Ferroelectrics -A Modern Perspective[END_REF] and it is impossible to resume all the conceptual steps which started since the definition of polarization.

One of the key points that one has to face in dealing with this delicate aspect is the common presence of PBC. The first problem that arises is the inherent discontinuity carried by the infinite replication of the simulation box. At the edges of the boxes infinite electric field strengths clearly occur when a linear electrostatic potential is applied within. However, this problem can be more or less easily handled by confining the MD with some specific constraint even if this solution is not so rewarding.

More seriously, the periodicity in the presence of a macroscopic electric field E will lead to a change in the electron potential in each replication of the simulation box. The intrinsic problem resides in the non-periodic nature of the position operator. In particular, the electric field changes by a factor of eE • R under a traslation by a lattice vector R and even a small field varies the nature of the energy eigenstates. Moreover, because the potential is unbounded from below, the ground state is ill-defined [START_REF] Nunes | Real-space Approach to Calculation of Electric Polarization and Dielectric Constants[END_REF][START_REF] Wannier | Wave Functions and Effective Hamiltonian for Bloch Electrons in an Electric Field[END_REF][START_REF] Nenciu | Dynamics of Band Electrons in Electric and Magnetic Fields: Rigorous Justification of the Effective Hamiltonian[END_REF].

Many perturbative treatments of the application of an electric field have been proposed but only within the Modern Theory of Polarization and Berry's phases this delicate issue can be efficiently tackled [START_REF] Resta | Chapter of Physics of Ferroelectrics -A Modern Perspective[END_REF][START_REF] Berry | Quantal Phase Factors Accompanying Adiabatic Changes[END_REF]. Few of them [START_REF] Nunes | Real-space Approach to Calculation of Electric Polarization and Dielectric Constants[END_REF] are founded on a Wannier-function-based solution to the finite-field problem that was not very useful in practice [START_REF] Resta | Chapter of Physics of Ferroelectrics -A Modern Perspective[END_REF]. By means of the use of Bloch's functions, Nunes and Gonze [START_REF] Nunes | Berry-phase Treatment of the Homogeneous Electric Field Perturbation in Insulators[END_REF] showed how the common perturbative treatments could be directly obtained from a variational principle based on minimizing an energy functional F of the following form

F = E KS ({ψ kn }) -E • P({ψ kn }) . (1.82)
E KS ({ψ kn }) represents the KS energy per unit volume as a function of all occupied Bloch's functions and P({ψ kn }) is the usual zero-field Berry's phase expression for the electronic polarization. In practice this equation has to be minimized with respect to all {ψ kn } in the presence of the field. This way also the Bloch's functions will become functions of E, and F develops into a function that implicitly depends on the field.

Within the KS scheme the Berry's phase polarization of the non-interacting KS system is in general not the correct one [START_REF] Gonze | Density-polarization Functional Theory of the Response of a Periodic Insulating Solid to an Electric Field[END_REF][START_REF] Gonze | Density-functional Theory of Polar Insulators[END_REF]. One way to circumvent this problem is to consider a generalized HK theorem in which n(r) and P(r) uniquely determine any ground state property [START_REF] Gonze | Density-polarization Functional Theory of the Response of a Periodic Insulating Solid to an Electric Field[END_REF][START_REF] Gonze | Density-functional Theory of Polar Insulators[END_REF].

Moreover, in the corresponding KS scheme one has to find an effective periodic potential and an effective electric field E KS that yield for a non-interacting system the same n(r) and P(r).

In this theory [START_REF] Gonze | Density-polarization Functional Theory of the Response of a Periodic Insulating Solid to an Electric Field[END_REF][START_REF] Gonze | Density-functional Theory of Polar Insulators[END_REF] the polarization is the correct one by construction, whereas a correction E KS -E has to be introduced (i.e., one has to introduce a sort of "exchange-correlation electric field").

Only in 2002, Umari and Pasquarello [START_REF] Umari | Ab initio Molecular Dynamics in a Finite Homogeneous Electric Field[END_REF] demonstrated that the functional (1.82) is suitable for use as an energy functional for a variational approach to the finite-field problem as well.

The proof is not obvious since the occupied wavefunctions {ψ kn } are not eigenstates of the Hamiltonian. Hence, the problem of computing the polarization P in an electric field provides from a certain point of view the solution to the problem of computing any property of an insulator in a finite homogeneous electric field; it is in fact the introduction of the Berry's phase [START_REF] Berry | Quantal Phase Factors Accompanying Adiabatic Changes[END_REF] polarization into the funcional F of eq. (1.82) that solves the problem. In particular, by considering an electric field along the x direction, the following variational energy functional can be builded

E E [{ψ i }] = E 0 [{ψ i }] -E • P [{ψ i }] , (1.83)
where E 0 [{ψ i }] is the well-known energy functional in the zero-field system and P [{ψ i }] is the polarization defined by Resta [START_REF] Resta | Quantum-mechanical Position Operator in Extended Systems[END_REF]:

P [{ψ i }] = - L π Im(ln detS[{ψ i }]) , (1.84)
where L is the periodicity of the cell and S[{ψ i }] is the following matrix

S i,j = ψ i |e 2πix/L |ψ i (1.85)
for the set of doubly occupied wavefunctions. By comparing their results with other perturbational approaches [START_REF] Nunes | Berry-phase Treatment of the Homogeneous Electric Field Perturbation in Insulators[END_REF], Umari and Pasquarello conclude that their method is applicable "even when the periodic symmetry is relaxed" [START_REF] Umari | Ab initio Molecular Dynamics in a Finite Homogeneous Electric Field[END_REF]. Moreover, this formulation allows not only to treat finite homogeneous electric fields in first-principles calculations but also can be extended to yield the atomic forces. Indeed by adding the following term to the functional (1.83)

E E ion =-E • P ion , P ion = N ion i=1 Z i • R i , ( 1.86)
where P ion is the ionic polarization, R i is the position coordinate in the field direction and Z i is the charge of the ionic core, this definition leads to an extra-term on the force acting on the ith atom equal to

F i = EZ i .
Finally, we want to remark that the briefly treated methodology developed by Umari and Pasquarello [START_REF] Umari | Ab initio Molecular Dynamics in a Finite Homogeneous Electric Field[END_REF] has been thoroughly employed in the simulations presented in this thesis when an electric field has been applied.

Collective variables and Metadynamics

As stated before, ab initio methods are by their nature an extremely powerful tool in predicting and analyzing chemical reactions. However, one of the most important parameters in describing activated processes is of course its free energy surface (FES). This latter is represented by a function that tipically depends on a number of variables (i.e., collective variables (CV)) which can be either few and relatively numerous depending on the specific reaction. There exist several theoretical and computational methods that are capable to reproduce the FES of a chemical process and few of them have been employed in this thesis.

Let us consider a system in the canonical ensemble; we introduce a collective variable q(R), a function of the atomic coordinates, that is able to distinguish the relevant metastable states of the system (i.e., reactants and products). The probability of finding the system in a specific configuration characterized by the reaction coordinate s is given by

P (s) = 1 Q e - U (R) k B T δ(s(R) -s)d N R . (1.87)
The free energy is related to this quantity by

F (s) = -k B T ln P (s) , (1.88)
and thus

P (s) ∝ e - F (s) k B T .
As a consequence, if it were be possible to explore the entire configuration space of a system, what is in principle feasible by means of an extremely long equilibrium Molecular Dynamics (MD) trajectory, it would be straightforward to reconstruct its FES. Unfortunately, the rates at which chemical reactions evolve are easily smaller than 1 hour -1 making infeasible the application of even classical MD techniques in most cases. However, enhanced sampling algorithms allow overcoming this limitation and reconstruct the FES: in this thesis we employed metadynamics (MetD) [START_REF] Laio | Escaping Free-Energy Minima[END_REF]. This approach biases the potential energy along a CV to enhance the sampling of, in principle, the corresponding free energy landscape. The bias is represented by a history dependent function which, as the dynamics evolve, decreases the probability of visiting configurations already explored. In particular, by depositing (tipically Gaussian) hills of potential energy centered in the visited points of the CV space, MetD is able to fill and escape the local free energy minima and, at the end of the simulation, to reconstruct the underlying FES. At time t, the bias at s is

W (s, t) = N hills (t) i=1 he -1 2 ( s-s(t i ) σ i ) 2 , t i = i • τ , (1.89)
where τ is the inverse of the hills deposition rate (i.e., the MetD time step).

Step by step the filling procedure of the "valleys" characterizing the FES will lead to a flattening of this latter once the convergence of the calculation has been achieved. Indeed a central assumption of MetD is that

-lim t→∞ W (s, t) ≃ F (s) + const. (1.90)
The statistical convergence of every MetD calculation is a key aspect in accurately reconstructing the correct FES. Although in expression 1.90 the concept of infinite limit has been introduced, for practical purposes this circumstance is achieved when the free energy profile show a parallel growth of the FES (flattening). This circumstance occurs when fixed-height Gaussian hills of potential are deposited but it exists also the possibility to employ the welltempered MetD variant [START_REF] Barducci | Well-Tempered Metadynamics: A Smoothy Converging and Tunable Free-Energy Method[END_REF]. In this approach, potential hills of progressively reduced heights are deposited when the simulation proceeds, becoming almost negligible once the convergence is achieved.

When MetD is exploited in conjuction with CPMD (due to the high computational cost of ab initio forces), the history dependent potential has to force the system to cross the energy barriers in few picoseconds. This implies that a huge amount of energy is injected in the sample in the form of bias potential and this may lead to a significant dishomogeneity in the temperature distribution and to instability in the dynamics. To address this problem an extended Lagrangian approach has been introduced [START_REF] Iannuzzi | Efficient Exploration of Reactive Potential Energy Surfaces Using Car-Parrinello Molecular Dynamics[END_REF], and auxiliary fictitious variables have been added to the CP Lagrangian:

L = L CP + Ns α=1 1 2 µ s α ṡ2 α (t) - Ns α=1 1 2 [S α (R(t)) -s α (t)] 2 + V (t, [s]) , (1.91)
where the second term is the fictitious kinetic energy associated with the fictitious variables {s α }. These variables are auxiliary degrees of freedom with associated fictitious masses µ s α which are introduced to explore the CV space. In this way, the dynamics of these extra-degrees of freedom can be explicitly controlled by suitable thermostats. An important requirements within CPMD is indeed the adiabatic separation with respect to the electronic degrees of freedom. If µ α is very large, the CV will relax to equilibrium distribution very slowly, reducing the efficiency of the method. Although its formulation appears to be elegant and conceptually generalizing, in practice the adiabatic separation can be maintained without the briefly presented extended Lagrangian approach, as it is in our cases [START_REF] Laio | Metadynamics: a Method to Simulate Rare Events and Reconstruct the Free Energy in Biophysics[END_REF].

As briefly stated in the Introduction, the necessity of the a priori choise of the CV which drive the reaction hugely limit the applicability of pure MetD techniques. Indeed it is not straightforward to predict which is the minimum set of CV that, describing the dynamical evolution of the system, are essential in reconstructing the FES. Moreover, the volume of the space of the CV grows exponentially with the number of CV making use of more than ≈ 2 computationally expensive.

Path Collective Variables with a new definition of distance

By the knowledge of the dynamical evolution of a given chemical reactioni.e., its trajectory -the problem of the a priori choise of a set of CV can, in principle, be bypassed. Branduardi et al. [START_REF] Branduardi; Gervasio | From A to B in Free Energy Space[END_REF] introduced the so-called "path-CV" which are constructed from a tentative reaction path {R 0 } connecting the reactants and the products states. By defining a measure D of the distance between the current reaction {R} and the reference one {R 0 } as, e.g., the mean square difference between all pairs of the interatomic distances computed in the reference states (say A and B for the reactants and the products states, respectively)

D({R A }, {R B }) = 2 N (N -1) N J>I (|R A I -R A J | -|R B I -R B J |) 2 , (1.92)
one can construct a CV that mark the progress of the reaction along the path:

s(t) = N f k=1 ke -λD(R(t),R k ) N f k ′ =1 e -λD(R(t),R k )
.

(

Of course, one can define the measure D in many different manners. Moreover, another important CV is introduced within this approach:

z(t) = - 1 λ ln N f k=1 e -λD(R(t),R k ) , (1.94)
which measures the distance away fromi.e., orthogonal to -the reference path.

The reference frame {R k } composing the tentative path can be built based on the available information about the system (i.e., transition state at T = 0), or by some interpolation between the end points. The big advantage of this approach is that the system is guided along the transition, still retaining the freedom to find a new path.

As an extension of path-CV, a new metric D suitable for chemical reactions has been very recently developed by Pietrucci and Saitta [START_REF] Pietrucci | Formamide Reaction Network in Gas Phase and Solution via a Unified Theoretical Approach: Toward a Reconciliation of Different Prebiotic Scenarios[END_REF].

These latter define a new metrical concept of the measure D which reads

D(R(t), R k ) = IS [C IS (t) -C k IS ] 2 , (1.95)
where C IS is the coordination number between atom I of species S ′ and all atoms J of species S whereas C k IS is its value referred to the reference structure {R k }. This definition of distance is well-suited to reactions since coordination numbers are permutation invariant with respect to atoms belonging to a given species S: any hydrogen, for instance, can take part in a protonation event (i.e., it can come from the solvent, from the solute, etc.). In practice the quantity C IS is defined by means of the switching function

C IS (t) = J 1 -R IJ (t) R 0 SS ′ N 1 -R IJ (t) R 0 SS ′ M , (1.96)
where R 0 SS ′ take into account the natural bond lengths between distinct atoms (see also Fig. 1.1). But the most striking advantage of this technique resides in the fact that in principle no insights might be known about the reaction under study. For instance, no predetermined trajectories are needed in order to produce a desired dynamical evolution of the atomic entities. Indeed a simple contact matrix description of the initial and of the end states is required and uniquely the structural formulas of these latter are mandatory. In practice, using the C k IS entries obtained from a short MD simulation of the end states can improve the resolution of the intermediate and transition states. Moreover, by its formulation, this new method represents a straightforward way to compute FESs in topologically equivalent CV spaces of gas phase and condensed phase chemical reactions. In this way, also a direct comparison between these two situations can be made highlighting the role of the solvent in assisting more or less complex reactions. This method has been employed in all the free energy calculations presented in this thesis.

In many cases MetD calculations converge and accurate evaluations of the free energies can be made directly [START_REF] Pietrucci | Formamide Reaction Network in Gas Phase and Solution via a Unified Theoretical Approach: Toward a Reconciliation of Different Prebiotic Scenarios[END_REF]. In more articulated reaction processes it has been proved that reaching the convergence is not trivial, neither using MetD with fixed Gaussian heights nor in the welltempered variant [START_REF] Barducci | Well-Tempered Metadynamics: A Smoothy Converging and Tunable Free-Energy Method[END_REF]. Despite this, "exploratory" simulations can be carried out in order to discover the position in the CV space of the relevant basins characterizing the FES and of the (putative) transition states connecting them. Once obtained the "map" of the free energy landscape it is possible to perform a series of parallel umbrella sampling (US) [START_REF] Torrie | Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling[END_REF] simulations in order to sample the reaction path and accurately evaluate the underlying FES. For the sake of completeness, a brief mention to the US technique follows.

Umbrella Sampling and the Weighted (and Dynamic) Histogram Analysis Methods

The sampling probability (1.87) can be modified by virtue of a biased potential added to the system which depends only on the CV. For the sake of simplicity now we will take only one CV but the following treatment can be straightforwardly extended to multi-dimensional FESs. The potential energy will be thus modified in U (R) → U (R) + V (s). Hence the biased distribution of the CV will be

P ′ (s) ∝ e - U (R)+V (s) k B T δ(s -s(R))d N R ∝ e - V (s) k B T P (s) , (1.97)
and the biased free energy F ′ (s) results to be

F ′ (s) = -k B T ln P ′ (s) = F (s) + V (s) + const. (1.98)
The additional constant is irrelevant in determing FESs since only energy differences are meaningful. By removing a posteriori the bias potential one can simply recover the unbiased FES F (s)

F (s) = F ′ (s) -V (s) + const. (1.99)
In practice, one of the most used bias potential is the harmonic one defined as

V (s) = 1 2 k(s -s 0 ) 2 , (1.100)
which depending on the strength of the constant k is able to sample even very difficult-to-sample free energy barriers. A biased MD with such a kind of bias potential is known as Umbrella Sampling (US). One has to keep in mind that for a two-dimensional casei.e., a FES that depends on two CV, s and z -the just defined harmonic potential has the form

V (s, z) = 1 2 k s (s -s 0 ) 2 + 1 2 k z (z -z 0 ) 2 .
Clearly in the one-dimensional case, the sampled distribution will be

P ′ (s) ∝ P (s)e - k(s-s 0 ) 2 2k B T , ( 1.101)
and the consequent biased free energy can be evaluated as (1.98), whereas by means of eq.

(1.99) one can recover the unbiased free energy.

Of course, in general when performing US simulations the bias produces an enhanced sampling in regions close to the minimum s 0 of the bias. In order to sample a relevant portion of the FES one has to combine multiple independent US simulations centered in distinct s i . Moreover, to the aim of gather enough statistical information, each sampled window has to overlap with the other nearest sampled regions as better as possible. This way, by applying the Weighted Histogram Analysis Method (WHAM) [START_REF] Kumar | The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method[END_REF][START_REF] Ferrenberg | Optimized Monte Carlo Data Analysis[END_REF][START_REF] Bennet | Efficient Estimation of Free Energy Differences from Monte Carlo Data[END_REF][START_REF] Shirts | Statistically Optimal Analysis of Samples from Multiple Equilibrium States[END_REF] it is possible to combine in a statistically well-founded manner the sampled windows in the CV space and to map the whole unbiased FES. In practice, if we have to collect N sims distinct US simulations, WHAM equations that are able to reconstruct the underlying free energy over the sampled (distinct) regions are

P (s) = N sims i=1 n i (s) N sims i=1 n i (s)e (F i (s)-V i,bias (s)) k B T , (1.102) F i (s) = -k B T ln S bins P (s)e - V i,bias (s) k B T . (1.103)
In these two equations, which have to be solved by iteration to self consistency, n i (s) is the number of counts in histogram bin associated with a given value of s. V i,bias is the biasing potential in the (US) simulation i and F i (s) is the free energy shift. Finally P (s), which is unknown as well as F i (s), is the best estimate of the unbiased probability distribution provided that a predetermined tolerance threshold has been imposed in the iterative procedure.

Although WHAM is the most used technique in order to combine multiple US simulations and thus obtain accurate free energy estimates, its reliability resides on the achievement of the equilibrium in each sampled window. Indeed each US simulation has to fully sample (in principle) the corresponding equilibrium probabilities. In practice this circumstance is not always verified leading to instabilities in the WHAM estimates of the FES [START_REF] Rosta | Free Energies from Dynamic Weighted Histogram Analysis Using Unbiased Markov State Model[END_REF]. Very recently, a dynamic version of this method (i.e., the DHAM) has been introduced [START_REF] Rosta | Free Energies from Dynamic Weighted Histogram Analysis Using Unbiased Markov State Model[END_REF]. DHAM has been proved to be more general than WHAM, with the latter representing only a special limiting case of its dynamical variant [START_REF] Rosta | Free Energies from Dynamic Weighted Histogram Analysis Using Unbiased Markov State Model[END_REF]. It makes use of a global Markov state model in order to obtain the free energy over a number of CV. In particular, a maximum likelihood estimate of the Markov transition matrix is constructed by joint unbiasing of the transition counts along discretized reaction coordinates. Finally, within DHAM an explicit approximation of the stationary distribution of the resulting Markov matrix has been introduced, which does not require the above presented iterative procedure characterizing WHAM.

In this thesis the US method has been used in order to converge the "exploratory" MetD sampling obtained by means of the new CV approach developed by Pietrucci and Saitta [START_REF] Pietrucci | Formamide Reaction Network in Gas Phase and Solution via a Unified Theoretical Approach: Toward a Reconciliation of Different Prebiotic Scenarios[END_REF].

At the end of this process, a recombination of the multiple US simulations has been performed via WHAM and thus the FESs reconstructed.

Chapter 2

Results on Proton Transfer in H-bonded Systems

Water ices under an electric field Introduction

Water is one of the most abundant substances on Earth and likely is the most studied compound in condensed matter physics. It has a rather baroque phase diagram with at least sixteen crystalline phases and three amorphous phases [START_REF] Petrenko | Physics of Ice[END_REF]. However, despite the huge experimental evidence and theoretical results which have accumulated so far, some structural and dynamical properties of the liquid and solid phases of water are not yet fully understood [START_REF] Petrenko | Physics of Ice[END_REF][START_REF] Debenedetti | Metastable Liquids -Concepts and Principles[END_REF].

Of paramount relevance is the phenomenon of molecular dissociation and the related phenomenon of proton transfer (PT). This process, which in liquid water is known as (auto)protolysis, plays a crucial role in many disparate domains, from neurobiology to electrolytic batteries and hydrogen-based technology [START_REF] Kaila | pH and Brain Function[END_REF][START_REF] Zoulias | Hydrogen-Based Autonomous Power Systems[END_REF]. At standard conditions this phenomenon is as rare as fundamental since, among other things, it represents the basis of the definition of the pH concept via the well-known reaction

2H 2 O ⇋ OH -+ H 3 O + , ( 2.1)
where two water molecules ionize by transferring one proton. Although the way in which this process -known as Grotthuss mechanism [START_REF] Grotthuss | Sur la décomposition de l'eau et des corps qu'elle tient en dissolution à l'aide de l'électricité galvanique[END_REF] -occurs along hydrogen-bonded chains of water molecules is phenomenologically clear since the '800s, a coherent theoretical framework which may explain in detail all the relevant aspects characterizing the protolysis and the PT mechanism is still missing.

In principle, the study of the H-bond network of solid water phases is easier than that of the liquid phase and hence ice phases are good candidates where to study the properties of the PT mechanism. At ambient conditions, the stable crystalline phase is ice I h , whose molecules are positionally arranged in an hexagonal lattice while being orientationally disordered (i.e., the proton sublattice is randomly distributed). Upon cooling KOH-doped ice I h down to a temperature of 72 K, it is possible to obtain its proton-ordered counterpart (ice XI) [START_REF] Tajima | Phase Transition in KOH-Doped Hexagonal Ice[END_REF][START_REF] Umemoto | Anomalous Pressure-Induced Transition(s) in Ice XI[END_REF],

whose molecules are oriented along the c-axis (i.e., the z-component of the dipole moments is everywhere positive). In the I h to XI first-order phase transition, the KOH hydroxide ions catalyze the rearrangement of H bonds, via defects formation, which leads to an enhanced molecular mobility. This way, water molecules exploit a more stable condition by reducing the original symmetry of the crystalline phase from hexagonal to orthorhombic, thus giving rise to the ferroelectric ice XI phase.

Because of the extremely rare nature of the protolysis phenomenon, Saitta and coworkers [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF] have shown that it is possible to stimulate and study the proton transport process in a systematic fashion by applying an external electric field.

In this Section, as already done in a correspective paper [START_REF] Cassone | Proton Conduction in Water Ices under an Electric Field[END_REF], we apply this computational approach to crystalline phases of water. We use the Car-Parrinello technique [START_REF] Car | Unified Approach for Molecular Dynamics and Density-Functional Theory[END_REF] to determine the energetic profiles of the molecular dissociation phenomenon and monitor some dynamic properties of the system, with particular emphasis on PT along H-bonded water molecules.

We also carry out some comparisons with the experimental results available in ices with high defects concentrations.

Finally, in order to shed more light on these apparently elaborate phenomena, we present a series of ab initio simulations that are performed to understand and exploit the symmetry differences between ice I h and ice XI. This aim is achieved by applying a static electric field along two different crystalline directions. Upon monitoring and analyzing the PT activity in the two topologically different materials, important conclusions can be drawn. This kind of analysis is of paramount importance because the hydrogen migration properties are relevant in many disparate fields [START_REF] Mcintosh | Quantum Effects in the Diffusion of Hydrogen on Ru[END_REF] and the related (an)isotropic piezoelectric effects are important [START_REF] Bystrov | Molecular Modeling of the Piezoelectric Effect in the Ferroelectric Polymer Poly(vinylidene fluoride) (PVDF)[END_REF] also for H-bonded systems [START_REF] Werling | Piezoelectric Hydrogen Bonding: Computational Screening for a Design Rationale[END_REF]. Moreover, here we investigate (as already done in a paper [START_REF] Cassone | Effect of Electric Field Orientation on the Mechanical and Electrical Properties of Water Ices: An Ab-initio Study[END_REF])

several dynamical properties of H-bonded molecular "wires" in the bulk during the proton migration process.

In the following Section the computational methods employed in these specific works concerning the ice phases are presented in detail.

Computational methods

By using the Car-Parrinello approach [START_REF] Car | Unified Approach for Molecular Dynamics and Density-Functional Theory[END_REF] we performed a series of static and dynamic ab initio simulations on two ice phases, i.e., ice I h and ice XI, under the effect of intense electric fields, applied along the c cell parameter, in order to study the phenomenon of ionic conduction. All data were obtained using the software suite Quantum ESPRESSO [START_REF] Giannozzi | QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulation of Materials[END_REF] for ab initio electronic structure calculations within which it is possible to handle macroscopic polarization and finite electric field by exploiting the modern theory of polarization (Berry's phases) [START_REF] Umari | Ab initio Molecular Dynamics in a Finite Homogeneous Electric Field[END_REF][START_REF] Berry | Quantal Phase Factors Accompanying Adiabatic Changes[END_REF].

Our systems were composed of 64 H 2 O molecules arranged in an orthorhombic cell with pa-rameters a = 9.05 Å, b = 15.67 Å, and c = 14.87 Å. In order to test potential size effects along the field direction, we also performed some numerical calculations on samples with 96 and 128 molecules, built with a c cell parameter whose values were 22.31 Å and 29.74 Å, respectively.

Periodic boundary conditions were applied in all simulations.

The positions occupied by the oxygen atoms are the same in the two ice structures; however, the arrangements of the hydrogen atoms differ significantly. In fact, the orientations of the molecules that compose the I h phase, while complying with the Bernal-Fowler ice rules [START_REF] Bernal | A Theory of Water and Ionic Solution, with Particular Reference to Hydrogen and Hydroxyl Ions[END_REF], are disordered, giving rise to a P6 3 /mmc symmetry; instead, those belonging to the ice XI phase are orientationally ordered, and form a crystalline structure that belongs to the Cmc2 1 space group.

The calculation of the molecular dissociation energy profiles is performed via electronic minimization by using, sequentially, the Gram-Schmidt and Davidson algorithms. In the following, we shall refer to such total energy profiles as "single", "double" or "triple" barriers according to the number of protons involved in the jump from their original sites. In the case of a double (triple) barrier, in which two (three) protons of H-bonded molecules simultaneously move, the protons were selected in this way: one belonging to the donor, one to the acceptor molecule, which in turn becomes the donor, and so on. We verified that the single and double-barrier cases cannot account for a potentially alternative scenario which may justify the jump of a proton on the energetic side. Hence, we limit our discussion to the triple-barrier profiles that are displayed in Fig. 2.1, whose inset also shows three snapshots of the molecular configurations involved in the calculations. These configurations were chosen because the molecular dipole moments are partially aligned with the field direction (positive z component). In particular, the configuration shown in Fig. 2.1-a is the same as the one used in the ferroelectric case and is obviously the most favored one for PT in the arrangement of the present numerical experiment. Barriers along the selected reaction coordinates were calculated via single, double or triple proton shift(s) of 0.038 Å.

The dynamic ab initio simulations were performed at the nominal temperature of 300 K which, according to the literature [START_REF] Del Buono | Model Dependence of Quantum Isotope Effects in Liquid Water[END_REF][START_REF] Allesch | First Principles Simulations of Rigid Water[END_REF][START_REF] Grossman | Towards an Assessment of the Accuracy of Density Functional Theory for First Principles Simulations of Water[END_REF], corresponds to a real temperature of about 250 K. The effect of the electric field was explored for intensities less or equal to 0.36 V/ Å, with a step increment of approximately 0.036 V/ Å , which leads to a total simulation time of more than ten picoseconds (for each value of the field intensity the duration of the simulation was at least 1 ps).

Several dynamic runs were conducted also for stronger electric fields (up to about 0.73 V/ Å), in order to obtain a full description of the current-voltage diagram of ice I h . The fictitious electronic mass was set at a value of 300 a.u., with a cutoff energy ranging between 25 and 35

Ry and a cutoff energy for the charge density of 200 -280 Ry: these values made it possible to adopt a Car-Parrinello timestep of about 0.1 fs. We used the Perdew-Burke-Ernzerhov (PBE) functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF], which belongs to the generalized gradient approximation (GGA), as an exchange and correlation functional because of its well known reliability when dealing with aqueous systems [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF][START_REF] Fiolhais | A Primer in Density Functional Theory[END_REF][START_REF] Santra | Hydrogen Bonds and van der Waals Forces in Ice at Ambient and High Pressure[END_REF][START_REF] Hamann | H 2 O Hydrogen Bonding in Density-Functional Theory[END_REF]. The electronic interaction with the nuclei was described through ultrasoft pseudopotentials (USPP) and the ion dynamics was carried out in the NVT ensemble, using the Verlet algorithm and Nosé-Hoover thermostat at a frequency of 13.5 THz. Protons were treated as classical particles. Although nuclear quantum effects do certainly play a role in proton barriers and diffusion, we focused more on relative values as emerge from a comparison between proton-ordered and proton-disordered phases.

The conductivities were obtained by exploiting Ohm's law. The current intensity is related to the number of charge carriers and to the time interval considered. The former quantity was obtained from the simulation data by calculating the number of protons travelling through a section area of the simulation box, orthogonal to the electric field direction. The other quantities, such as the potential difference in the simulation box and the correspoding resistance/conductivity, can be computed in a straightforward way.

Once the general properties of the PT mechanism in the abovementioned ice phases have been elucidated, a detailed investigation of the mechanical and the electrical responses is proposed under the action of intense electric fields applied along two directions, viz., the c-axis and the b-axis. For each value of the field, we run the dynamics for about 2 ps, cumulating a total simulation time longer than 20 ps for each sample and field direction. Although all the calculations performed here exploit the PBE functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF], in the light of some recent debate on its reliability for aqueous systems [START_REF] Santra | Hydrogen Bonds and van der Waals Forces in Ice at Ambient and High Pressure[END_REF][START_REF] Morales | Quantum Monte Carlo Benchmark of Exchange-Correlation Functionals for Bulk Water[END_REF][START_REF] Alfè | Analyzing the Errors of DFT Approximations for Compressed Water Systems[END_REF], we carry out some calculations with the BLYP functional [START_REF] Becke | Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior[END_REF] as well. The energies of the distorted structures and the corresponding equations of state (EOS) are also calculated using the same DFT-based Quantum ESPRESSO code. Following the energy-strain approach, we execute a series of relaxation simulations by using the steepest-descent technique combined with the damping dynamics of the ions for different values of the cell parameters, while keeping the cell volume fixed over the whole range of explored field strengths. With this approach it is also possible to calculate, inter alia, the c 66 component of the stiffness tensor (i.e., one of the elastic constants), that is related to the curvature of the EOS around its minimum. In hexagonal systems five stiffness constants only are independent, viz., c 11 , c 12 , c 13 , c 33 , and c 44 , since c 66 = 1 2 (c 11 -c 12 ) [START_REF] Landau | Theory of Elasticity[END_REF]. It is possible to determine this latter quantity through the relation (expressed here in the Voigt matrix notation)

c 66 = 1 V ∂ 2 E ∂ε 6 ∂ε 6       E=E 0 , (2.2)
where V is the volume, ε denotes an infinitesimal deformation, and E is the self-consistent total energy of the system. In practice, the derivative is carried out on a quadratic fit of the data for the energy.

General results and discussion

In liquid water molecules start to dissociate for field intensities of about 0.25 V/ Å [START_REF] Geissler | Autoionization in Liquid Water[END_REF][START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF][START_REF] Stuve | Ionization of Water in Interfacial Electric Fields: An Electrochemical Overview[END_REF][START_REF] Hammadi | Proton and light ion nanobeams from field ionization of water[END_REF],

while a net proton flow is observed along the z-direction for values of the intensity above 0.36 V/ Å [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF]. We found that the same thresholds apply to ice I h which, in this respect, behaves as liquid water; this finding is not particularly surprising since the dipole moments orientational distributions are similar in both phases. Instead, in ice XI the angular symmetry is broken in that the z-component of the dipole moments is, on average, positive, a condition that leads to the emergence of an internal electric field (ferroelectric phase). As a result, the electrostatic coupling with the external field is significantly enhanced. Correspondingly, the first dissociation events are observed for field intensities in the range 0.18 V/ Å -0.22 V/ Å, at which threshold protons also start to flow in a sustained way. The simultaneous onset of such two phenomena is a distinguishing feature of ice XI as compared with the two-stage behavior observed in ice I h and liquid water.

In order to investigate the mechanism underlying molecular dissociation and PT we first performed a series of "static" calculations of the potential-energy profiles, obtained in both ice structures at 0 K after displacing -along the field direction and by the same distance from their original sites -up to three protons, while keeping the oxygen atoms fixed at their original positions: one selected proton belonging to a molecule which possesses at least one hydrogen atom which can be shifted along the positive direction of the field, a proton covalently bonded to a first-neighbor oxygen atom, and one more proton belonging to a molecule that is hydrogen-bonded to the second one. In ice XI all choices are equivalent since they at most differ for the orientation of the central molecule with respect to the other two; instead, in ice I h there are three configurations only of three "wired" molecules such that on each molecule it is possible to find at least one hydrogen atom, bonding a neighboring oxygen, which can be displaced along the field direction (i.e., along the positive z axis). The results of these static calculations are shown in Fig. 2.1 for a field strength of 0.36 V/ Å, which corresponds to the proton conduction regime. We see that, in addition to the energy minimum associated with the unperturbed configuration, all the curves show a secondary minimum at the distance where the displaced proton might covalently bond with the nearby oxygen atom. However, even for very high field strengths this minimum is not deep enough to justify, on a mere energetic basis, the PT and the corresponding formation of an ionic pair. Hence, a scenario where oxygen atoms are held fixed while protons move along a "chain" does not suffice to explain the dissociation of a water molecule. Note, however, that, notwithstanding the failure of our constrained model calculation in accounting for the energetics of the PT, the height of the barrier separating the two minima in the most favorable I h configuration (∼ 1.4 eV) turns out to be in fair agreement with the available experimental data for the activation energy of the process leading to the formation of an ionic pair [START_REF] Petrenko | Physics of Ice[END_REF].

The role of the network formed by the oxygen atoms in assisting, locally, the PT event can be elucidated by dynamical simulations. We observed that the external field produces a notable effect in that it "flattens" the planes parallel to its direction. This phenomenon is more clearly resolved in ice XI, as shown in the inset of Fig. 2.2, and can be assimilated to a sort of "pressure" [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF][START_REF] Schwegler | Dissociation of Water Under Pressure[END_REF] exerted by the applied field which forces the mutual alignment of molecular dipole moments. The resulting polarization effect ultimately enhances the electrostatic coupling between the molecular dipole moments and the external field in a powered manner. In ice I h the flattening of the planes was not observed as clearly as in the ferroelectric case, because a larger energy is needed, on average, to align disordered dipole moments with the field direction. In ice XI we further observed the alignement of the oxygen planes perpendicular to the applied field: molecular orientations are restored after a diffusion process has taken place via a Grotthuss mechanism [START_REF] Grotthuss | Sur la décomposition de l'eau et des corps qu'elle tient en dissolution à l'aide de l'électricité galvanique[END_REF], which causes an inversion of the involved dipole moments.

As a consequence of the flattening process described above, the mutual distances between oxygen atoms decrease: this structural rearrangement enhances the probability of occurrence of correlated proton jumps which enable the phenomenon of diffusion. Correspondingly, the potential barrier that protons must overcome when migrating from a donor to an acceptor molecule, following the formation of a Zundel complex (H 3 O + +H 2 O→H 5 O 2 + ), decreases [START_REF] Marx | Proton Transfer 200 Years after von Grotthuss: Insights from Ab initio Simulations[END_REF] and attains its minimum value for an oxygen-oxygen distance of about 2.4 Å [START_REF] Xie | The Protonated Water Dimer: Extensive Theoretical Studies of H 5 O + 2[END_REF][START_REF] Tuckerman | On the Quantum Nature of the Shared Proton in Hydrogen Bonds[END_REF]. In both crystalline structures we observed a shortening of the hydrogen-bond chain whenever an ionization event occurred (see Fig. 2.2). This behavior is coherent with the phenomenology observed by Hassanali and coworkers [START_REF] Hassanali | On the Recombination of Hydronium and Hydroxide Ions in Water[END_REF] in liquid water and confirms that the diffusion event is the result of a highly cooperative molecular motion. When the

H 3 O + -H 2 O distance
decreases, the proton localization probability decreases as well and, correspondingly, the PT probability increases along the hydrogen-bond chain [START_REF] Marx | Proton Transfer 200 Years after von Grotthuss: Insights from Ab initio Simulations[END_REF]. The process is Zundel-to-Zundel type and is mediated by water molecules as in the liquid phase [START_REF] Hassanali | On the Recombination of Hydronium and Hydroxide Ions in Water[END_REF]: when an H 3 O + ion gets sufficiently close to an H 2 O molecule, a Zundel ion forms and, after multiple exchanges of the ionic roles between the two entities, the complex eventually propagates with a Grotthuss-type dynamics.

As soon as a sustained proton diffusion has been activated, both ices display an ionic ohmic behavior (see Fig. 2.3), which has been observed experimentally in the common phase of ice for voltages up to about 10 V [START_REF] Maidique | Transfer of protons through pure ice in single crystals -Part III. The dielectric relaxation spectra of water, ice, and aqueous solutions and their interpretation[END_REF]. We calculated the static conductivities of ice I h and ice XI which turned out to be 6.9 Ω -1 cm -1 and 5.8 Ω -1 cm -1 , respectively. Indeed, these values are many orders of magnitude higher than those obtained experimentally for pure crystalline phases of ice which, however, refer to significantly lower field strengths (10 -10 V/ Å-10 -4 V/ Å) [START_REF] Maidique | Transfer of protons through pure ice in single crystals -Part III. The dielectric relaxation spectra of water, ice, and aqueous solutions and their interpretation[END_REF][START_REF] Liu | Determination of Phase Diagram of Water and Investigation of Electrical Transport Properties of Ices VI and VII[END_REF]. Following Jaccard who explained proton mobility in ices in terms of defects [START_REF] Jaccard | Thermodynamics of Irreversible Processes Applied to Ice[END_REF], we attribute the discrepancy between our estimates and the experimental data to the largely different concentrations of ionic defects that are generated by the applied field. Moreover, as shown in Fig. 2.3, no significant and systematic differences emerge between the estimates obtained with different sample sizes. Table 2.1 shows the percent fractions of defects that are present in the two simulated crystalline structures for increasing field strengths; the values corresponding to the most intense applied field (0.73 V/ Å) range between 16% and 19%. Such concentrations are actually typical of doped materials. Grimm et al. [START_REF] Grimm | Low-Frequency Electrical Properties of Polycrystalline Saline Ice and Salt Hydrates[END_REF] investigated ice I h doped with CaCl 2 and other salts for increasing molarities (1 M, 2.25 M, 3.14 M, and 3.96 M), over a wide range of temperatures.

In fact, the concentration of defects generated via the formation of HCl (see Table 2.2) with a salt molarity of 3.96 is comparable with the values displayed in Table 2.1. The estimates of the percent fractions of defects present in experimentally heavily-doped ices were obtained by considering a one-to-one correspondence between doping and defect levels, as shown in [START_REF] Fletcher | The chemical physics of ice[END_REF].

Upon extrapolating the low-temperature conductivity data of this material outside the sampled rangei.e., well above the eutectic point, up to the temperature of our numerical experiments -one would recover values of the same order or even larger than those reported above for non-doped ice.

Another comparison can be made with the conductivity of KOH-doped ice I h as measured by Zaretskii et al. [START_REF] Zaretskii | The Protonic Conductivity of Heavily KOH-doped ice[END_REF] with alternate currents. However, beyond a temperature threshold of about 200 K such a frequency-dependent conductivity does no longer depend, in an appreciable way, on the frequency of the applied field. For instance, in 1 M KOH-doped ice I h , where the fraction of ionic defects is about 2%, the conductivity measured by Zaretskii et al. at 283 K was about 10 -1 Ω -1 cm -1 . Again, the discrepancy between this value and our results for the static conductivities of non-doped ices can be largely accounted for by the different temperatures and defects concentrations. Indeed, upon extrapolating the experimental conductivity curves, we recover values whose orders of magnitude are comparable with the data extracted by the present simulations.

Electric field strength (V/ Å)

Mechanical properties

In c). However, the most significant difference is found in the curvature: in ice XI the shape of the EOS changes with the field strength when the field is oriented along the c-axis, a behavior that, instead, is not exhibited by ice I h as well as by ice XI under the effect of a field parallel to the b-axis. whose projection onto the field direction varies from negative to positive values, prevent a softening of the two structures since a conspicuous molecular re-orientation of the dipoles is impossible in the mentioned crystalline phases during the relaxation experiment. We finally note that the dielectric constants of the two samples, when calculated with the PBE and with other exchange and correlation functionals [START_REF] Casazza | Ab initio Investigation of Electronic and Vibrational Contributions to Linear and Nonlinear Dielectric Properties of Ice[END_REF][START_REF] Murray | Dispersion Interactions and Vibrational Effects in Ice as a Function of Pressure: A First Principles Study[END_REF][START_REF] Lu | Dielectric Properties of Ice and Liquid Water from First-Principles Calculations[END_REF], exhibit very similar values; hence, we may argue that they do not influence in an appreciable way the ambivalent behavior of the mechanical response discussed above.

Electrical properties

We had already presented a series of results about the proton static conductivities of ices I h and XI under the action of an external electric field oriented along the c-axis [001] direction [START_REF] Cassone | Proton Conduction in Water Ices under an Electric Field[END_REF].

We had interpreted the differences in the respective dissociation and current thresholds by considering the different orientational distributions of the water molecules in the two crystalline structures. In addition, we had explained the resulting static conductivities through the percent fractions of field-induced defects. Here we want to discuss the intrinsic electrical anisotropy that is due to the proton arrangements by orienting the electric field along another direction of the crystal, i.e., the b-axis [010] direction, in the common hexagonal phase of ice as well as in its proton-ordered counterpart. The close correspondence of the dissociation (0.25 V/ Å) and current (0.36 V/ Å) thresholds in ice I h and liquid water [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF][START_REF] Stuve | Ionization of Water in Interfacial Electric Fields: An Electrochemical Overview[END_REF] was explained through the similarity of the distributions of dipole moments in the two phases.

However, in the proton-ordered case we found substantially similar threshold values of the field strength for both processes of molecular dissociation and sustained proton diffusion (0.18 -0.22 V/ Å). We noted that these latter estimates fall in a range of values lower than those found in ice I h and liquid water as a consequence of the enhanced coupling between the external field and the dipole moments; moreover, since the polarization effects are more important in the ferroelectric case, the local electrical environment enhances the probability of molecular dissociation and protonic diffusion along the field direction thus giving rise to a single-threshold process, at variance with the I h and liquid water cases. The PT mechanism occurs along a sequence of H-bonds until a proton has travelled a certain path which becomes infertile as for the transfer of another proton. The re-activation of a given proton transmitting H-bonds passage is achieved via the formation of Bjerrum defects [START_REF] Petrenko | Physics of Ice[END_REF][START_REF] Jaccard | Thermodynamics of Irreversible Processes Applied to Ice[END_REF].

Upon orienting the external electric field along the b-axis direction, it is possible to compare and understand the coupling between dipole moments and the external electric field, thus obtaining partial insights on the time-dependent local polarization effects. On account of the quasi-random molecular orientation distribution in common ice, one would expect a cartesian axis reversal invariance with respect to a change in the field propagation direction. However, in the ferroelectric phase, a different coupling -as compared with the case shortly described above (electric field along the ferroelectric axis) -arises from the intrinsic electrical anisotropy of the material [START_REF] Tajima | Phase Transition in KOH-Doped Hexagonal Ice[END_REF][START_REF] Umemoto | Anomalous Pressure-Induced Transition(s) in Ice XI[END_REF][START_REF] Jackson | Thermally-Stimulated Depolarization Studies of the Ice XIâĹŠIce Ih Phase Transition[END_REF]. events occur for E=0.22 V/ Å (the process is signalled by the absolute minimum pointed by the arrow); b) average molecular distances in an H-bond wire in ice XI under the action of a field oriented along the c-axis and in which the first dissociation events occur for E=0.18 V/ Å; the intermolecular length shows a huge increment after attaining its absolute minimum.

is found in the molecular dissociation and sustained proton diffusion threshold values, for any orientation of the field; instead, in the ice XI samples these thresholds are modified in a sensitive way and a two-stage process takes place.

In ordinary ice the ionization event becomes manifest for 0.22 V/ Å, but a net proton flow occurs for field intensities of 0.32-0.36 V/ Å. We note that, notwithstanding the fact that a net proton flow can be measured for a field strength of about 0.32 V/ Å, the discrepancy with the result obtained after orienting the field along the c-axis falls within the error of our calculations; moreover, this discrepancy can be related to the small enhancement of the dielectric anisotropy introduced by the DFT exchange and correlation functional used [START_REF] Schönherr | Dielectric Properties of Water Ice, the Ice I h /XI Phase Transition, and an Assessment of Density Functional Theory[END_REF].

The field strength threshold of the dissociation events in ice XI is also registered, somewhat surprisingly, for 0.22 V/ Å , but the most favourable circumstance after molecular dissociation is the rattling of the proton shared between the two forming ions, and only for an intensity of 0.29 V/ Å we do observe the first net proton flow events. This two-stage process is the first neat alteration of the response behavior as compared to the case in which the field direction is parallel to the ferroelectric axis. We also note that the corresponding molecular dissociation thresholds in the simulations carried out with the BLYP functional are 0.25 V/ Å and 0.22 -0.25 V/ Å for ice I h and ice XI, respectively. Differences with respect to the PBE results fall within the error range of our calculations, and the estimates of protonic current thresholds are the same for the two functionals.

Dynamics of H-bond chains

With regard to the local dynamics, when the dissociation events occur we observe a shortening of the H-bond lengths that allows for a fast PT along the respective reaction coordinates, as shown in Fig. 2.6 for the two crystalline phases. In the ice XI case the average oxygen-oxygen distances are slightly shorter than in the hexagonal common phase, which suggests a more efficient transfer process of defects.

We define the H-bond wire length as the sum of the individual donor-acceptor distances between H-bond connected oxygen atoms in which a molecular dissociation has occurred. Figure 2.7-a shows the behavior of an ice I h H-bond wire, whose length has been rescaled by dividing it by the number of involved intermolecular distances and in which a combined PT between the oxygens has occurred in a few hundreds femtoseconds, followed by a ready proton recombination. A pronounced relative minimum is visible just before the global one indicating a first (unsuccessful) attempt at molecular dissociation and local PT. We observe another curious aspect which characterizes all the investigated H-bond wires: after dissociation, the length of the "wire" always increases with respect to the average one and attaines a global maximum just after the global minimum, showing a sort of chain "breath"; in all cases (even in the ferroelectric one-stage mechanism of PT shown in Fig. 2.7-b) subsequent correlated molecular dissociations, produced immediately later along the same H-bond path, do not occur. This finding is consistent with the Jaccard theory of defects propagation [START_REF] Jaccard | Thermodynamics of Irreversible Processes Applied to Ice[END_REF]. Hence, we argue that the mechanism underlying molecular dissociation, which implies a certain degree of cooperation between the oxygen atoms, is the same in liquid water [START_REF] Hassanali | Proton Transfer Through the Water Gossamer[END_REF][START_REF] Hassanali | On the Recombination of Hydronium and Hydroxide Ions in Water[END_REF] as well as in the crystalline phases that we have investigated.

Ice I h vs. ice XI: anisotropy effects

As already noted, the field intensity value that allows for a sustained protonic current in common ice is comparable to that observed in liquid water (0.36 V/ Å) [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF][START_REF] Cassone | Proton Conduction in Water Ices under an Electric Field[END_REF][START_REF] Stuve | Ionization of Water in Interfacial Electric Fields: An Electrochemical Overview[END_REF] and does not substantially depend on the field direction. The response of ice I h to static electrical perturbations is intrinsically isotropic on account of the topological isotropy conferred to this ice phase by its proton sublattice. Incidentally, this result recalls the hypothesized isotropic nature of the G correlation parameter (also known as dipolar autocorrelation coefficient) in ice I h ; this finding is also supported by Monte Carlo simulations and is intimately related to the electrical response via the dipole moments distribution [START_REF] Petrenko | Physics of Ice[END_REF].

However, when the field is oriented along the b-axis of ice XI, the behavior of the material changes significantly if compared with that of the same structure under an electric field oriented along the c-axis. In these conditions the molecular dissociation events and the subsequent proton diffusion are not simultaneous, as previously noticed for ice XI under the action of a field parallel to the ferroelectric direction. In particular, as the concentration of both types of defects becomes relevant, the system starts to conduct in almost the same way as in the case of ice XI when the field lies along the ferroelectric axis, evolving, from a certain point of view, to a partial (unstable and field-induced) ferroelectricity.

In order to give a more complete description of the ionic conduction in these systems, we also calculate the ionic current-voltage diagrams corresponding to our unit cells. The results are displayed in Fig. 2.8-a for the two phases when the electric field is oriented along the b-axis; c). This result can be interpreted at relatively low voltages (up to about 5 V for our unit cell) by considering that the system under the action of a field oriented along the b-axis, generates a sufficient concentration of defects, up to a field strength of 0.29 V/ Å (ionic current threshold). Above this value, the sample is somewhat indistinguishable, as for the electrical response, from a similar sample under the effect of an electric field acting along the ferroelectric axis. The hypothetic slight differences in the derivatives (i.e., in the conductivities)

Partial ferroelectricity along the b-axis of ice XI

When a sample of ice XI is subject to a field strength equal to or greater than 0.29 V/ Å along the b-axis, it is possible to observe a sort of partially ferroelectric domains along the field axis in the presence of a conspicuous amount of defects. The process can be explained as follows: starting from the symmetry of the dipole moments distribution of the ice XI molecular structure, one can observe that there are some planes (as the middle plane almost parallel to the xy plane in the orthorhombic centered unit cell shown in Fig. 2.9) in which molecules Also ice I h has a certain number of molecules that are mostly oriented along the b-axis with a similar y component of their dipole moment as reported above, but there are only triplets of connected molecules with this peculiarity, whose chain extremities are H-bonded to other chains of H-bonded molecules that consistently exhibit lower (and in most cases negative) ycomponents of their dipole moments, and so on. In this way, when a field is applied along any direction, the surrounding environment prevents any potential ferroelectric transition in the common ice phase, even for high field intensities.

These observations become crucial when a static electric field is applied along the b-axis. In fact, the first molecules which start to dissociate in ice XI are those belonging to the planes mentioned above and, as the concentration of defects increases with the external field intensity, the system develops a partial ferroelectricity along the Hamiltonian discontinuity axis along which the field is oriented. This behavior is likely due to the collective positive dipole moment (with respect to the field direction) exhibited by these planes, which is enhanced by the field itself. This leads, in turn, to an enhancement of the local electric field which produces more and more defects and, correspondingly, favors a larger molecular mobility. These circumstances are accompanied by a shortening of the transition time [START_REF] Tajima | Phase Transition in KOH-Doped Hexagonal Ice[END_REF]. It is also known that the H-bond network is slightly distorted by defects [START_REF] Petrenko | Physics of Ice[END_REF]. The H-bonds that lie transversal to the field direction get weaker and weaker as the field intensity increases [START_REF] Shevnukov | Electric field induced transitions in water clusters[END_REF]. All these considerations justify the fact that the molecules have the opportunity of slightly rotating around the more robust H-bond, making a statistical attempt at the formation of Bjerrum defects. The creation of ionic propagating defects, beyond the current threshold, is not by itself sufficient to explain the ionic conductivity. The formation of Bjerrum defects is crucial in that they can "open" new protonic paths along H-bonds [START_REF] Petrenko | Physics of Ice[END_REF][START_REF] Jaccard | Thermodynamics of Irreversible Processes Applied to Ice[END_REF] after an ion has diffused away. The propagation of defects takes place via the Grotthuss mechanism [START_REF] Grotthuss | Sur la décomposition de l'eau et des corps qu'elle tient en dissolution à l'aide de l'électricité galvanique[END_REF].

Summing up, different pieces of evidence -such as the basic PT mechanism, the field-induced enhancement of molecular mobility, the generation of ionic defects, and the stronger local polarization effects created in one half of the ab planes -make the formation of ferroelectric patches possible in ice XI under the effect of an electric field with a strength larger than 0.29 V/ Å and projected along the b-axis.

This phase is obviously different from the standard ice XI structure as one can see from the distribution of H-bonds (i.e., the simmetry along the b-axis is different from the simmetry along the c-axis). However, by looking at the problem simply from the electrical point of view, the two conductors are practically the same. In fact, also in the ice XI case subject to a field along the c-axis and above the current threshold, it is possible to find molecular configurations in which molecules are obviously present that have a negative z-component of the dipole moment just after/before being involved in PT: this behavior is inherent to the conduction process generated by protonic defects.

It is important to remark that a complete ferroelectric phase transition is impossible; in fact, in such a case, the system would melt because of the formation of alternating series of D-and L-defects between those planes that we have formerly called ferroelectric and anti-ferroelectric (with respect to the b-axis) (see Fig. This finding is not particularly striking if one imagines that also in a ferroelectric sample under the effect of a field parallel to the c-axis at some instants some molecules are oriented oppositely with respect to the field in order to eventually "open" some H-bond paths for another PT or to "close" them after a PT has taken place along a given path.

Liquid methanol under an electric field Introduction

Methanol (CH 3 OH) is the simplest aliphatic alcohol. Single-donor H-bonds, confer upon this system an important role for the understanding of the proton transfer (PT) mechanism which can be basically schematized according to the following reaction:

2CH 3 OH ⇋ CH 3 O -+ CH 3 OH 2 + . (2.3)
Indeed, because of the nature of the methanol molecule, its H-bond organization can be seen as a set of cluster structures characterized by different topologies [START_REF] Pauling | The Nature of the Chemical Bond[END_REF][START_REF] Sarkar | Molecular Clusters and Correlations in Liquid Methanol at Room Temperature[END_REF][START_REF] Montague | Structural Studies of Liquid Alcohols by Neutron Diffraction[END_REF][START_REF] Narten | Hydrogen Bonding in Liquid Methanol[END_REF][START_REF] Magini | On the Structure of Methyl Alcohol at Room Temperature[END_REF][START_REF] Wilson | X-ray Absorption Spectroscopy of Liquid Methanol Microjets: Bulk Electronic Structure and Hydrogen Bonding Network[END_REF][START_REF] Jorgensen | Quantum and Statistical Mechanical Studies of Liquids. 11. Transferable Intermolecular Potential Functions. Applications to Liquid Methanol Including Internal Rotation[END_REF][START_REF] Haughney | Molecular-dynamics Simulation of Liquid Methanol[END_REF][129][START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF][START_REF] Morrone | Ab initio Molecular Dynamics Study of Proton Mobility in Liquid Methanol[END_REF].

Liquid methanol displays a fast PT and incidentally this process takes place along the H-bond structure via the Grotthuss mechanism [START_REF] Grotthuss | Sur la décomposition de l'eau et des corps qu'elle tient en dissolution à l'aide de l'électricité galvanique[END_REF]. In this latter one assumes that a protonic defect migrates along an H-bond network through a series of correlated hoppings. As a consequence, the orientation of the molecular dipole vector changes soon after a PT has occurred. A deeper understanding of the PT phenomenon is of crucial importance since it represents the basis of the operational mode of many biologically and technologically important systems [START_REF] Kaila | pH and Brain Function[END_REF][START_REF] Zoulias | Hydrogen-Based Autonomous Power Systems[END_REF].

The H-bond structure of liquid methanol has been extensively studied experimentally [START_REF] Pauling | The Nature of the Chemical Bond[END_REF][START_REF] Sarkar | Molecular Clusters and Correlations in Liquid Methanol at Room Temperature[END_REF][START_REF] Narten | Hydrogen Bonding in Liquid Methanol[END_REF][START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Weitkamp | Hydrogen Bonding in Liquid Methanol at Ambient Conditions and at High Pressure[END_REF][START_REF] Tsuchida | Density-functional Study of Liquid Methanol[END_REF][START_REF] Torrie | Structure of the α-phase of Solid Methanol[END_REF][START_REF] Tauer | On the Crystal Structures, Residual Entropy and Dielectric Anomaly of Methanol[END_REF] as well as with theoretical [START_REF] Pauling | The Nature of the Chemical Bond[END_REF][START_REF] Suresh | Influence of Electric Field on the Hydrogen Bond Network of Methanol[END_REF] and computational methods [START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF][START_REF] Boyd | A Density Functional Study of Methanol Clusters[END_REF][START_REF] Wang | A Hybrid Quantum Mechanical Force Field Molecular Dynamics Simulation of Liquid Methanol: Vibrational Frequency Shifts as a Probe of the Quantum Mechanical/Molecular Mechanical Coupling[END_REF].

Tetrahedrally coordinated networks are present in water whereas methanol molecules form linear chains which are interrupted by bifurcations [START_REF] Narten | Hydrogen Bonding in Liquid Methanol[END_REF][START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Weitkamp | Hydrogen Bonding in Liquid Methanol at Ambient Conditions and at High Pressure[END_REF]; an enormous effort has been done in investigating methanol clusters [START_REF] Boyd | A Density Functional Study of Methanol Clusters[END_REF]141].

The H-bond network certainly plays a crucial role in assisting the PT phenomenon both in water [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF][START_REF] Agmon | Elementary Steps in Excited-State Proton Transfer[END_REF] as well as in methanol [START_REF] Morrone | Ab initio Molecular Dynamics Study of Proton Mobility in Liquid Methanol[END_REF][START_REF] Chang | Isomeric Transitions between Linear and Cyclic H + (CH 3 OH) 4,5 : Implications for Proton Migration in Liquid Methanol[END_REF] but the microscopic comprehension of this phenomenon is not fully clear yet. To obtain a more realistic picture of the process which underlies the PT activity, several attempts have been made by simulating, with ab initio methods, pure liquid methanol with an excess proton [START_REF] Morrone | Ab initio Molecular Dynamics Study of Proton Mobility in Liquid Methanol[END_REF] and methanol-water solutions [START_REF] Morrone | Ab initio Molecular Dynamics Simulation of the Structure and Proton Transport Dynamics of Methanol-Water Solutions[END_REF]. The response to an external electric field and the consequential rearrangements of methanol clusters have also been investigated [START_REF] Bounds | A Molecular Orbital Study of the Effect of an External Electric Field on Methanol[END_REF][START_REF] Rai | Methanol Clusters (CH 3 OH) n , n = 3 -6 in External Electric Field: Density Functional Theory Approach[END_REF].

Methanol is used in apparatus for the production of electrical energy based on the chemical reaction cited above. One example is provided by direct methanol fuel cell (DMFC) technology where electrical energy is produced by transferring protons from the alcoholic group of methanol molecules through a membrane [START_REF] Dohle | [END_REF]. Hence, a deeper understanding of the microscopic mechanism that leads to the molecular dissociation and/or to a net proton flow is of paramount importance to increase the efficiency of several ecological systems which produce energy from methanol.

In the following, as already presented in a recent paper [START_REF] Cassone | Liquid Methanol under a Static Electric Field[END_REF], we present the first, to the best of our knowledge, ab initio molecular dynamics study of bulk liquid methanol under the effect of a static electric field. By applying an external field it is possible to study, in a systematic way, the PT phenomenon and to disclose the relevant reaction coordinates that underlie molecular dissociation and ionic (protonic) conduction processes. Analogies and differences with the corresponding phenomena in liquid water are also highlighted.

Theory and simulation

We used the software package Quantum ESPRESSO [START_REF] Giannozzi | QUANTUM ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulation of Materials[END_REF], based on the Car-Parrinello (CP) approach [START_REF] Car | Unified Approach for Molecular Dynamics and Density-Functional Theory[END_REF], to perform ab initio calculations of a sample of methanol molecules under the action of intense electric fields applied along a given direction (corresponding to the z axis).

Our sample contained 32 CH 3 OH molecules (i.e., 192 atoms) arranged in a cubic cell with side parameter a = 12.93 Å , so as to reproduce the experimental density of 0.79 g/cm 3 at about 300 K. As is customary, the structure was replicated in space by using periodic boundary conditions.

We carried out ab initio simulations at the nominal temperature of 300 K. We gradually increased the intensity of the electric field with a step increment of about 0.051 V/ Å from zero up to a maximum of 0.51 V/ Å. In the zero-field case we executed a dynamics of almost 6 ps whereas, for each other value of the field intensity, we ran the dynamics for about 2 ps, thus cumulating a total simulation time longer than 25 ps.

The fictitious electronic mass was set to a value of 300 a.u., with a cutoff energy of 35 Rydberg (Ry) and a cutoff energy for the charge density of 280 Ry, which allowed us to adopt a timestep of 0.12 fs. With such cutoff values the sample is described in a realistic way since the core electronic interaction is being depicted through ultrasoft pseudopotentials (USPP). As for exchange and correlation effects, we adopted the Perdew-Burke-Ernzerhof (PBE) functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF],

which belongs to the generalized gradient approximation (GGA) regime. This choice was motivated not only by the adherence of previous results [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF] with the available experimental data for systems under a static electric field [START_REF] Stuve | Ionization of Water in Interfacial Electric Fields: An Electrochemical Overview[END_REF][START_REF] Hammadi | Proton and light ion nanobeams from field ionization of water[END_REF][START_REF] Rothfuss | The Influence of the Surface Electric Field on Water Ionization: a Step Dissociative Ionization and Desorption Mechanism for Water Ion Cluster Emission from a Platinum Film Emitter Tip[END_REF] but also by the good performance of the PBE functional in H-bonded systems because of an appropriate description of polarization [START_REF] Marom | Dispersion Interactions with Density-functional theory: Benchmarking Semiempirical and Interatomic Pairwise Corrected Density Functionals[END_REF]. The dynamics of ions was simulated classically within a constant number, volume, and temperature (NVT) ensemble, using the Verlet algorithm and a Nosé-Hoover thermostat set at a frequency of 13.5 THz.

We performed a structural analysis of the H-bond network and identified an H-bond through the following geometric conditions (that should be simultaneously fulfilled) [START_REF] Haughney | Molecular-dynamics Simulation of Liquid Methanol[END_REF][START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF]

: (i) r(O• • • H) bond distance smaller than 2.6 Å; (ii) O• • • O distance smaller than 3.5 Å; (iii) ∠HO• • • O bond
angle smaller than 30 • . We calculated also the time autocorrelation function of H-bonds as:

C I (t) = i,j s ij (t 0 )s ij (t 0 + t) i,j s ij (t 0 ) , (2.4)
where the indices i and j run on all pairs of first-neighbor molecules which at t 0 were Hbonded, t 0 being the time at which the measurement process begins; s ij = 1 if the criterion for the presence of an H-bond is fulfilled, s ij = 0 otherwise. The results were averaged over different initial configurations.

We finally computed the conductivity σ using Ohm's law in the field regime where this law is satisfied by our model liquid.

Results

Average instantaneous correlations and H-bond dynamics

The In the past decades the local structure of pure liquid methanol has been studied by neutron diffraction [START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Adya | The Structure of Liquid Methanol by H/D Substitution Technique of Neutron Diffraction[END_REF], X-ray scattering [START_REF] Narten | Hydrogen Bonding in Liquid Methanol[END_REF], and DFT-level simulations [START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF][START_REF] Morrone | Ab initio Molecular Dynamics Study of Proton Mobility in Liquid Methanol[END_REF][START_REF] Tsuchida | Density-functional Study of Liquid Methanol[END_REF][START_REF] Handgraaf | Density-functional Theory-based Molecular Simulation Study of Liquid Methanol[END_REF].

The accuracy of the ab initio methods in reproducing the average molecular arrangements has been confirmed by several works, such as those by Morrone et al. [START_REF] Morrone | Ab initio Molecular Dynamics Study of Proton Mobility in Liquid Methanol[END_REF] and Pagliai et al. [START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF], who used the Becke-Lee-Yang-Parr (BLYP) functional. In this work we used instead the PBE functional, whose reliability has been already tested for water [START_REF] Fiolhais | A Primer in Density Functional Theory[END_REF][START_REF] Santra | Hydrogen Bonds and van der Waals Forces in Ice at Ambient and High Pressure[END_REF][START_REF] Hamann | H 2 O Hydrogen Bonding in Density-Functional Theory[END_REF]; the results listed in Table 2.3 for the position of the first peak of atomic RDFs show that this functional is a very good choice for liquid methanol as well; indeed, the present PBE results compare with the experimental data more favorably than the corresponding BLYP estimates [START_REF] Morrone | Ab initio Molecular Dynamics Study of Proton Mobility in Liquid Methanol[END_REF] or the estimates obtained with Monte Carlo simulations [START_REF] Shilov | Long Range Order and Hydrogen Bonding in Liquid Methanol: A Monte Carlo Simulation[END_REF]. Moreover, the O-O RDF determined by Pagliai et al. [START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF], which is in excellent agreement with the experimental data of Yamaguchi et al. [START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF], was obtained for a liquid density of 0.888 g/cm 3 , significantly higher than the experimental value (0.79 g/cm 3 ) at 298 K. Hence, we surmise that the less enhanced liquid structure observed by Morrone et al. is not likely to occur due to the relatively small value used for the truncation of the plane-wave series (25 Ry) but more plausibly is an artifact of the BLYP functional. Rather surprisingly, Handgraaf et al. [START_REF] Handgraaf | Density-functional Theory-based Molecular Simulation Study of Liquid Methanol[END_REF] recovered a good agreement of some RDFs with neutron diffraction data [START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF] by using the same setup of Pagliai et al. [START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF] (i.e., BLYP functional, Martins-Troulliers pseudopotentials, and a plane-wave basis cutoff of 70 Ry), but for a liquid density of 0.791 g/cm 3 .

Another characteristic of interest that can be obtained from the O-H RDF is the mean number of H-bonds per molecule n HB , which was evaluated by integrating the RDF up to the position of the first minimum. The results are in fair agreement with those obtained with the procedure suggested by Rapaport [START_REF] Rapaport | Hydrogen Bonds in Water: Network Organization and Lifetimes[END_REF], in which one separately determines the numbers of methanol molecules that are involved in zero-H, single-H, double-H, and triple-H-bonds. By adopting the geometric criterion for the identification of an H-bond in liquid methanol [START_REF] Haughney | Molecular-dynamics Simulation of Liquid Methanol[END_REF][START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF], the calculation leads to the following H-bond percent fraction distribution in the absence of an electric field: f 0 = 6.25%, f 1 = 12.5%, f 2 = 68.75%, and f 3 = 12.5%, which imply a mean value of H-bonds n HB = 1.88, in excellent agreement with other two independent numerical estimates obtained, respectively, for deuterated systems (1.89) [START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF] and with a classical Molecular Dynamics study (1.902) [START_REF] Sun | Molecular Dynamics Simulation of Liquid Methanol under the Influence of an External Electric Field[END_REF] as well as with the available experimental data (1.9) [START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF].

The reliability of the PBE functional in the absence of an external field is a good premise for testing this functional and the changes in the average local structure of the system it accounts for when a static and uniform electric field has been switched on. We first run a zero-field CP dynamics about 5.5 ps long at 300 K and then applied a field of increasing intensity. As shown in Fig. 2.14, pair correlations between oxygen atoms are initially enforced; both maxima and minima of the O-O RDF are enhanced and sharpened; in other words, the liquid gets more ordered at short and medium distances. However, as the field intensity grows, this trend changes in a significant way: in fact, the height of the first peak, after reaching a maximum for E ≈ 0.15 V/ Å (see the inset of Fig. 2.14), exhibits a rather abrupt drop of about 50% for field intensities between 0.25 V/ Å and 0.35 V/ Å, and stays approximately constant thereon.

Correspondingly, the first dip rises as also does the second peak. The first peak and the immediately following dip of the O-H RDF behave in a similar way (see Fig. 2.15). In addition, the first coordination shell moves to moderately shorter distances and a weak shoulder gradually emerges at about 1 Å, which indicates an increasing freedom of hydrogens in sampling shorter relative distances from an oxygen atom. A further consequence of the application of the electric field on the system is the gradual aligning of increasing fractions of molecular dipole moments along the field direction. Let θ be the angle formed by dipoles with the z axis; as can be seen from Fig. 2.17, in the absence of the field the distribution P (θ) of molecular dipole orientations exhibits a maximum for θ = 90 • [START_REF] Suresh | Influence of Electric Field on the Hydrogen Bond Network of Methanol[END_REF].

However, when an electric field has been switched on, the maximum of the distribution shifts to lower values of θ since molecules tend to align with the field axis, the larger the number the stronger the field.

The perturbing effect of too strong electric fields on the stability of the H-bond network can be further appreciated through the behavior of the time autocorrelation function defined by Eq. 2 and plotted in Fig. 2.18. In the absence of an external field, the function decays rather slowly and, notwithstanding the time limitations of CP simulations, we can infer an average lifetime of H-bonds of the order of 10 ps or even longer, whereas for a field strength of 0.36 V/ Å this value dramatically drops to 2 -3 ps.

Summing up, the structural and dynamical evidence presented above clearly indicate that the effects of a stationary and uniform electric field on the stability of the H-bonds formed by methanol molecules in the liquid phase critically depend on the intensity of the field: moderately high fields enforce and stabilize the H-bond network while too strong fields eventually break the bonds and disrupt the network. This phenomenology can be explained by considering that molecular dipoles tend to align with the field direction, a circumstance which enhances the formation probability of H-bonds [START_REF] Suresh | Influence of Electric Field on the Hydrogen Bond Network of Methanol[END_REF]. However, for very high field intensities a full alignment of molecular dipoles with the field direction is not possible because of the topological constraint imposed by the location of the H-bond donor/acceptor sites on the oxygen atoms. The resulting geometrical frustration likely weakens the structural stability of the H-bond network.

Molecular dissociation and protonic current

The average dipole moment of methanol molecules in the liquid phase has been already calculated with ab initio methods implementing the BLYP functional [START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF]. Depending on the number of H-bonds in which a methanol molecule is involved, its dipole moment can assume values equal to 2.06 D (no H-bonds), 2.24 D (one H-bond), 2.71 D (two H-bonds), and 2.99 D

(three H-bonds), corresponding to an average value of 2.64 D [START_REF] Pagliai | Hydrogen Bond Dynamics in Liquid Methanol[END_REF]. Our result for this latter quantity is 2.65 D. A similar calculation for water molecules leads to an average value of 2.95 D [START_REF] Boero | Hydrogen Bonding and Dipole Moment of Water at Supercritical Conditions: A First-Principles Molecular Dynamics Study[END_REF][START_REF] Boero | Water at Supercritical Conditions: A First-Principles Study[END_REF][START_REF] Silvestrelli | Water Molecule Dipole Moment in the Gas and in the Liquid Phase[END_REF][START_REF] Silvestrelli | Electronic, and Bonding Properties of Liquid Water from First Principles[END_REF], higher than for methanol molecules. Correspondingly, the molecular dissociation threshold of liquid methanol turns out to be higher than that of liquid water [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF] as well as of water ices [START_REF] Cassone | Proton Conduction in Water Ices under an Electric Field[END_REF][START_REF] Cassone | Effect of Electric Field Orientation on the Mechanical and Electrical Properties of Water Ices: An Ab-initio Study[END_REF]. In fact, although we observed an increased vibration amplitude of the hydroxyl group for (and above) a field strength of about 0.26 V/ Å, indicating that some molecules are preparing for ionization, the first dissociation events were recorded for field strengths of 0.31 V/ Å and beyond. It has been shown that in liquid and solid water the oxygen atoms play a crucial role in facilitating the PT events by reducing their mutual distances [START_REF] Cassone | Proton Conduction in Water Ices under an Electric Field[END_REF][START_REF] Hassanali | On the Recombination of Hydronium and Hydroxide Ions in Water[END_REF]. In methanol a direct analysis of the O-O distances revealed that the oxygen atoms involved in a PT event spend, on average, much more time closer to their respective partners in order to transfer the proton. This means that the ionic complex The increase of the rate of H-bond breaking and reforming processes (i.e., of fluctuations in the H-bond network) due to the field cannot lead in methanol to a PT enhancement comparable with that of water. In fact, the average coordination number of a neutral methanol molecule is 2 and the same value is obviously found for the protonated species CH 3 OH 2 + . Hence, larger H-bond fluctuations do not give rise to a significant increase of the PT activity as in liquid water. This is the main reason for the flow of a feeble protonic current even for intense field strengths (see Fig. 2.20). In fact, although the onset of a net proton flow has been observed for the same field intensity recorded for liquid water (≈ 0.35 V/ Å) [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF], the efficiency of the process is extremely lower than in water: very weak currents of the order of tenths of µA are induced by fields with a strength of 0.36 -0.51 V/ Å. The resulting conductivity has an average value of 0.40 S cm -1 , which is more than one order of magnitude lower than that exhibited by liquid and solid water under the action of fields with similar intensities [START_REF] Saitta | Ab initio Molecular Dynamics Study of Dissociation of Water Under an Electric Field[END_REF][START_REF] Cassone | Proton Conduction in Water Ices under an Electric Field[END_REF][START_REF] Cassone | Effect of Electric Field Orientation on the Mechanical and Electrical Properties of Water Ices: An Ab-initio Study[END_REF]. Ionic conductivities of such an order of magnitude are involved in DMFC devices [START_REF] Duankaew | Methanol Permeability and Proton Conductivity of Direct Methanol Fuels Cell Membranes based on Sulfonated Poly (Vinyl Alcohol)-Layered Silicate Nanocomposites[END_REF] and are, on average, one order of magnitude higher than those tipically involved in modified Nafion membranes [START_REF] Kabumoto | Acta of the Conference of the Electrochemical Society on Energy/Battery Technology[END_REF].

[C 2 H 6 O 2 H 3 ] + (see

Chapter 3

Results on Methanol Chemical Reactions

Methanol chemistry at high electric fields

We report on an ab initio Molecular Dynamics study of bulk liquid methanol under the effect of strong static electric fields. We found that an unknown chemical reaction occurs leading to the production of formaldehyde and methane. Moreover, other reactions are observed by further increasing the field strength yielding, inter alia, dimethyl ether. Incidentally the two microscopic reaction mechanisms that give rise to formaldehyde and methane on the one hand and to dimethyl ether on the other, involve the presence of the methanol cation CH 3 OH 2 + and methoxide CH 3 O -as precursor species. To the aim to clarify the contributions of the atomic reactivities in assisting the electrostatic effects during this double channel reaction, a Fukui analysis and a Noncovalent Interactions study have been conducted on the two ionic species, also with respect to their neutral parent. In particular we found an appreciable increase of the carbon reactivity when considered in the cationic structure which, under some (local) circumstances, may lead to the release of CH 3 + as a transient species. Depending on the topological molecular arrangement of the methoxide in the neighborhood of the CH 3 OH 2 + cation, the reaction can evolve by following two drastically different paths: one leading to the formaldehyde and methane formation and the other producing dimethyl ether.

Field-induced chemical reactions A brief overview

Althuogh the electrochemistry of methanol should be an almost definitively explored topic, some chemical reactions which occur when an external electric field is applied on the bulk liquid are still unknown. At relatively weak field strengths it is possible to induce the cleavage of the alcoholic covalent bond which results in the concept of pH of methanol. Recently, it has been possible to quantify some of the parameters involved in this process; indeed, by applying on a sample of pure liquid methanol a static electric field equal to or greater than 0.30 V/ Å, it is possible to induce molecular dissociations in which the alcoholic proton can migrate via a Grotthuss-like mechanism [START_REF] Grotthuss | Sur la décomposition de l'eau et des corps qu'elle tient en dissolution à l'aide de l'électricité galvanique[END_REF] along the H-bond network [START_REF] Cassone | Liquid Methanol under a Static Electric Field[END_REF]. Of course this process is assisted by a certain fraction of ionic species, such as the methanol cation CH 3 OH 2 + and the methoxide CH 3 O -, which transiently exist and are responsible for the (ionic) charge flow in the system.

But, what happens at higher field intensities? Interestingly, at field intensities stronger than 0.50 V/ Å it is possible to observe several field-induced chemical reactions. As a consequence, at these field strengths the carbon-carbon radial distribution function exhibits the onset of a new first peak at shorter distances than the standard one, as shown in Fig. 3.1. Moreover, from this latter it is also possible to notice that the peak becomes more and more pronounced by further increasing the field strength. The appearance of a new typical length in the system is the signature of the birth of different chemical species stemming from methanol.

A field strength of 0.55 V/ Å, assisted by not negligible local (i.e., solvent) advantageous circum- stances, is able to break some covalent bonds in the methyl groups of few methoxide molecules and to create formaldehyde molecules. The formation of this species is of course accompained by the release of water and methane according to the following reaction

2CH 3 OH E -→ H 2 CO + CH 4 + H 2 O . ( 3.1)
To our best knowledge this is an unknown chemical reaction which involves three of the most widely employed molecular species in industries, such as methanol, methane, and formaldehyde. This latter is an extremely reactive compound that gives rise to a progressive complexification of the species in the system. Going in order of increasing complexity, another observed molecule is a consequence of the formaldehyde hydration:

H 2 CO + H 2 O E -→ CH 2 (OH) 2 , ( 3.2)
which is known as methanediol, formaldehyde monohydrate, or methylene glycol.

By experiencing the field action, methanol molecules can also formally react to each other in order to give rise to dimethyl ether species (also known as methoxymethane, CH 3 OCH 3 ) accordingly to the reaction

2CH 3 OH E -→ CH 3 OCH 3 + H 2 O, (3.3)
which leads to the release of water molecules. Incidentally, as it will explained later, this kind of reaction can also proceed through the intervention of the already present formaldehyde molecules but the most likely circumstance is so far the peculiar recombination process which takes place between the intermediate species CH 3 OH 2 + and CH 3 O -.

Finally, the biggest molecule produced by the action of the field is the so called methoxymethanol, which is known also as methanolic formaldehyde or formaldehyde methyl hemiacetal.

Indeed formaldehyde, as an extremely reactive compound, can directly react with the sorrounding methanol molecules. The nominal reaction which yields methoxymethanol can be expressed as

H 2 CO + CH 3 OH E -→ CH 3 OCH 2 OH. (3.4)

Analysis of the chemical reactions mechanisms

As already explained [START_REF] Cassone | Liquid Methanol under a Static Electric Field[END_REF], at low/mild-field regime (i.e., ≈ 0.30 V/ Å) despite obvious polarization effects, it is possible to observe a series of proton transfers. Indeed, intensities of this magnitude are able to break the alcoholic covalent bond and release protons in the sample.

Hence, above this threshold, the field effects will be manifest also as a consequence of the presence of several ions. This leads to locally enhanced contributions to the electrostatic potential which in turn produces an increase of the molecular reactivities.

Here we want to outline several reaction mechanisms that are observable at our thermodynamic conditions, leaving the detailed (atomic) reactivity study to the homonym section. For the sake of clarity we treat here the chemical reactions in a no-chronological way, introducing them in order of increasing importance of the respective reaction products and not in order of occurrence during the numerical experiment which, in turn, was conducted by increasing step by step the field intensity. The reader should only to keep in mind that the first nontrivial chemical reaction (i.e., the dissociation of the alcoholic hydrogen can be regarded as a trivial reaction) that occurs in the sample is the formation of formaldehyde, methane, and water, observable at a field strength of 0.55 V/ Å. Hence, starting from this threshold the original sample of pure methanol is also characterized by the presence of these species.

Methoxymethanol

One of the simplest reaction mechanisms leads to the formation of the most complex species observed in our sample: the methoxymethanol molecule. The first appearance of this latter occurs at a field strength of 0.60 V/ Å, when are already present molecules such as water, formaldehyde, and methane, besides the obvious methanol and its cationic and anionic counter- Although the interesting behavior of the species involved in this process, which react all together in concert, the reaction product is not a molecule of particular scientific or industrial interest.

H H O H H C O H H C O H H H H C O H E field direction H H O H H C O H H C O H H H H C O H H H O H H C O H H C O H H H H C O H Methoxymethanol
It is in fact a transient equilibrium species found usually in mixtures of formaldehyde and methanol, which may also contain water. This compound is not produced as a pure chemical entity, since it is normally unstable, and its primary application is as chemical intermediate in the manifacture of various urea formaldehyde and melamine formaldehyde resins [161].

Dimethyl ether

The formation of a species of great relevance has been observed in our numerical experiment: the dimethyl ether. It represents the simplest ether and it is usually employed as a chemical reagent or as a propellant for aerosol. Even more interestingly, it is used as "clean" fuel because of the reduced amounts of NO x and CO as products [162]. Commercially, methoxymethane (one of the several names of dimethyl ether) is produced via dehydration of methanol in the presence of a compound which acts as a catalysti.e., alumina - [START_REF] Tokay | Dimethyl Ether Synthesis over Alumina Based Catalysts[END_REF][START_REF] Xu | Synthesis of Dimethyle Ether (DME) from Methanol over Solid-Acid Catalysts[END_REF] whereas in our numerical experiment the catalyst of the chemical reactions is obviously the electric field.

There exist two main processes which lead to the onset of dimethyl ether; one of these, involving also a formaldehyde, can be depicted and dissected in few steps as shown in Fig. 3 formaldehyde and a methanol cation (left), the positive charge migration toward the field direction (center) leads to the formation of a dimethyl ether cation with the release of water (right). b) By obtaining an hydride from a methanol molecule (left), a dimethyl ether is formed (center) whereas another formaldehyde is created via the deprotonation of the ionic methanol species (center and right).

hyde along the field direction (Fig. 3.3-a, left), the positive charge transfer is energetically less favorable if it occurs via a cation rotation in order to donate its proton to the aldehyde. In fact, since the redistribution of the electronic density is a faster process than a whole molecular re-orientation and because of the energetic contribution carried by the electric field at this strengthi.e., 0.60 V/ Åthe positive charge can migrate along the cation itself reaching the methyl group. Moreover, this process is enhanced by the presence of the nearby formaldehyde which, being highly polarized and oriented by the field action, attracts the positive charge in a not negligible way. The concerted action of the aldehyde and of the external field, together with the surrounding presence of other ionic species, leads to the cleavage of the CO bond of the methanol cation (Fig. 3.3-a, center). As a consequence, a methenium CH 3 + is released and, after the well-known umbrella inversion, approaches the partially negative head of the formaldehyde molecule, i.e., the oxygen atom, leading to the formation of a dimethyl ether cation with It should be clear that the local environment plays a crucial role in assisting all the presented chemical reactions since it acts as a reservoir of proton accepting/donating sites. Moreover, when the field strength is intense enough, some hydrides, always stemming from methyl groups, act in a hyper-reactive manner in the system. It is important also to remark that the above discussed "transition states" have been in truth hyperdissected and the concerted motions shown in Fig. 3.3 occur in a hundred of femtoseconds.

E field direction H H O C H H H O C H H H H H O C H H H O C H H H Umbrella inversion H H O O C H H H C H H H Dimethyl ether
The just described methodology by which a dimethyl ether molecule rises is not the only one and another reaction channel, even simpler and more favored, can be followed in the system. Indeed, as a first consequence of the presence of an electric field, some molecules/ions tends to aligne their less electronegative sites toward the field direction as it is the case for the methoxide depicted in Fig. 3.4 (left). On the other hand, because of the activated proton transfer some methanol cations can be instantaneously found with their positively charged part against the field direction. As a consequence, at 0.60 V/ Å it is possible to observe methanol cations and anions neighbour pairs (Fig. 3.4, left) oriented each other in such a way that, with respect to the field direction, different alternating slices (orthogonal to the field) of positively and negatively charged portions of space can be detected. The neutralization process can hence take place as depicted in Fig. 3.4. The positive charge carried by the proton belonging to the methanol cation can reach the methyl group of its own molecule (Fig. 3.4, left). In this configuration the combined action of the external field and of the neighbour negative head of the methoxide cleaves the CO bond of the cation and results in the release of a methenium (Fig. 3.4, center). After the renowned "umbrella inversion", this transient species approaches the negatively charged methanol ion and a new species of dimethyl ether (with the release of water) can be detected in the sample (Fig. 3.4, right).

Again, this kind of process occurs in an handful of decades of femtoseconds. Since it is clearly simpler than the mechanism described in Fig. 3.3 and because of the nature of the sample, is the most likely reaction channel for the onset of dimethyl ether.

Formaldehyde and methane

Among the species formed in our sample, one of the most important is formaldehyde prominently but not only because of its hypotesized key role in the prebiotic world [START_REF] Cleaves | The Prebiotic Geochemistry of Formaldehyde[END_REF]. It is the simplest aldehyde and an intermediate oxydation state one carbon molecule. The relevance of this molecule lies in its potential prebiotic role suggested by its enormous chemical reactivity.

Formaldehyde readily undergoes a variety of addition and redox reactions to give products of biological importance such as sugars (i.e., in the formose reaction) and amino acids [START_REF] Cleaves | The Prebiotic Geochemistry of Formaldehyde[END_REF]. In particular, although there exists a not negligible amount of uncertainty on the efficiency of the formose reaction [START_REF] Shapiro | Prebiotic Ribose Synthesis: a Critical Analysis[END_REF], it is still believed that H 2 CO may be a precursor to ribose and other sugars, thus having given rise to the RNA or pre-RNA world.

In despite of the prebiotic implications which can rise from the fact that formaldehyde is at the same time the smallest aldehyde and an incredibly reactive species, this molecule covers a wide commercial importance in the world industrial economy. In this respect, the relevance of formaldehyde is indisputable since it has been estimated that only in the U.S. and in Canada over 4 millions of workers are directly or indirectly involved in its business. The products of H 2 CO industry are pervasive and they generate a substantial volume of sales, provide a sizable number of jobs, and contribute to the local economies in countless visible and not-so-visible ways. Moreover, despite for the Second World War period, the economy related to formaldehyde has steadly grown with time [START_REF]The Economic Benefits of Formaldehyde to the United States and Canadian Economies[END_REF]. Industrially, the conversion of methanol in formaldehyde covered a big slice of this business but the known transformation processes employed were operative at gas-phase conditions [START_REF] Hader | Formaldehyde from Methanol[END_REF]. For all these reasons, a deep comprehension of the ways by which this fundamental molecule can originate is hence of primary interest.

Incidentally, the otherwise unknown chemical reaction that we have found (reaction 3.1) involves also another product of unquestionable industrial interest: the methane. Although the usual practice is to reduce methanol to methane, a single (i.e., one step) chemical reaction involving these two species with formaldehyde has never been observed, nor in nature, nor in a laboratory.

Here we want to show the basics of this reaction that is ubiquitous in our sample and it is also the most favorable one; it occurs in fact for the first time at the lowest of the field intensities able to induce some chemical variations (i.e., 0.55 V/ Å).

Because of the field-induced ionizations of the methanol molecules -at and above a field intensity of 0.30 V/ Å -the system is characterized also by the presence of two ionic species: CH and CH 3 O -, i.e., the methanol cation and the methoxide, respectively. Although the strength of 0.55 V/ Å is surely considerable, the role of the solvent appears to be central in the observed typical reaction mechanism for the birth of formaldehyde and methane species stemming from liquid methanol. Indeed in a condensed phase, at a given instant, are present several molecules arranged in such a way that the most electronegative "pole" is oriented in favor of the field direction (e.g., see the methoxide in the left side of Fig. 3.5). Moreover, also as a consequence of the activated proton transfer which occurs along the H-bond network, several methanol cation may be oriented with the excess proton in opposition to the electric field direction (Fig. 3.5, left). This circumstance, as already briefly explained, reaches its minimum energy configuration through the cleavage of the CO bond and the release of a methenium as a transient species.

Such a process is observable because this recombination process involves a re-arrangement of the electronic structures rather than a full re-orientation of both molecules which would led to a proton transfer. Of course, the methenium release is strongly encouraged by the mechanical stress exerted on the CO bond of the CH . This species will be employed for the formation of a dimethyl ether molecule (b). These sketches show also the key role played by the solvent in assisting these chemical reactions as a reservoir of proton and hydride accepting/donating sites (b). The overall process occurs in about 300 fs.

The transition state is characterized by the known "umbrella inversion" of the methenium just after the rupture of the CO bond and the release of an hydride stemming from the methyl group of the methoxide (Fig. 3.5, center). Hence, the recombination process between these two transient species leads to the formation of a methane whereas a formaldehyde and a water molecule appear all together as the final and stable equilibrium state (Fig. 3.5, right).

The chemical reaction just qulitatively described represents a new link between three important species: methanol, formaldehyde, and methane. This latter molecule is substantially inert and it does not participate in any other chemical reaction. On the other hand, the oxidation of methanol to formaldehyde leads to an increase of the molecular reactivity in the system. Indeed, the newly created aldehyde will be employed in the formation of many more complex species under the action of the field, as also described in the previous two Subsections.

What is really important is that the mechanism through which formaldehyde forms in our numerical experiment is in practice that one just analyzed. The reported example of the release of formaldehyde via the interaction between a dimethyl ether cation and a methanol molecule (Fig. 3.3-b) is statistically irrelevant with respect to the mechanism depicted in Fig 3 .5. This is likely a consequence of the nature of our original system: a sample of pure liquid methanol.

In order to show the multiplicity and the field-induced cooperation between the species, it is also possible to observe a mixture of the processes described in Fig. 3.3 and 3.5, as shown in From the left side of Fig. 3.6-a it is clear that the formaldehyde (and methane) formation via the ionic channel interaction picture does not depend on the relative orientation of the methanol cation and anion with respect to the field orientation, since in Fig. shown in the left side of Fig. 3.4. On the other hand, if the closer parts of the two neighbouring ions are represented by the respective methyl groups, the reaction will evolve in the formation of formaldehyde, methane, and water. This double channel reaction is schematically depicted in Fig. 3.7. These latter two important and still unknown reactions will be the topic of the next Sections.

In particular, since the transient ionic species of methanol appear to play a central role in the synthesis and in the catalysis of several crucial species, the focus will be pointed on their reactivity at an atomic scale which will be also compared with that of the neutral methanol molecule.

Stability of the formed species

At this point it is relevant to clarify the stability of the species in the system. The numerical experiment has been conducted by progressively increasing the field strength with a step increment of 0.05 V/ Å from zero up to a maximum of 0.75 V/ Å. In the zero-field case a dynamics of almost 6 ps has been ran whereas, for each other value of the field intensity, a dynamics of at least 2 ps has been performed, thus cumulating a total simulation time of almost 40 ps.

Moreover, in order to test the stability of the created species also in absence of the field, a dynamics longer than 10 ps for each just formed molecular species in the sample has been conducted by turning off the external field.

At the end of the investigation performed when the field intensity has been rised up to 0.55 V/ Å, in the sample was present the equivalent of ≈ 3 % of water, ≈ 3 % of methane, and ≈ 3 % of formaldehyde; this circumstance corresponds in having ≈ 6 % of dissociated methanol molecules if compared with the number of methanol molecules which were present at the beginning of our numerical experiment. Therefore, the percent fractions shown in Table 3.1, are determined as the ratio between the number of molecules of a given species and the number of total methanol molecules composing the sample at the beginning of the simulation (i.e., 32 molecules=192 atoms).

When the field intensity has been increased up to 0.60 V/ Å and 0.65 V/ Å, the formaldehyde molecules created in the meanwhile have been employed in the formation of dimethyl ether, formaldehyde monohydrate, and methoxymethanol, yielding the percent fractions of these species shown in Table 3.1 (second data column). These values correspond to have had ≈ 22 % of methanol molecules which have been formally involved in the abovementioned chemical reactions via molecular dissociation. For higher field strengths (i.e., 0.70 and 0.75 V/ Å) the system is characterized by an elevated degree of mixture of species since it is possible to observe the presence of all of the species consequence of a precursor water molecules release.

As briefly abovementioned, another kind of study on the molecular stabilities has been conducted. Indeed by turning off the external field, at a field strength of 0.55 V/ Å, once the first observed formaldehyde, water, and methane molecules have been formed in the system, they rest totally stable. In fact, by performing a dynamics longer than 10 ps, during which the (un)protonated methanol species neutralize by recombining, the species of water, methane, and formaldehyde rest totally stable. Only the water molecule has been involved in some of the proton transfer processes necessary for the neutralization of the ionic species.

Finally, it must be stressed that by switching off the electric field when all the studied neutral species have been already appeared in the systemi.e., above 0.60 V/ Åall of them rest stable, suggesting that the free energy minima in which they have fallen is deep enough at our thermodynamic condition (i.e., T=300 K).

Reactivity and (regio)selectivity

Since the most important chemical reactions observed in our numerical experiment are those leading to the formation of H variations in the global and regional reactivities of these ions with respect to their neutral parent will be treated.

However, before to do that, a preliminary analysis of the intermolecular interactions between the local reactants within the Noncovalent Interactions (NCI) framework could be of interest.

Indeed, as briefly abovementioned, the two reactions proceed also without any external application of the electric field in the gas phase and the path followed in each process is strongly dependent on the initial configuration of the two ions, as shown in Fig. 3.7. By considering the reduced density gradient after a single point calculation performed via the PBE [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] exchange and correlation (XC) functional, the NCI have been evaluated, as shown in Fig. 3.8.

This kind of calculation has been conducted by following the standard procedure [START_REF] Contreras-García | NCIPLOT: a Program for Plotting Noncovalent Interaction Regions[END_REF] for the evaluation of the reduced density gradient and therefore by selecting a correct isosurface in order to map the real space region in which the NCI acts. As it is clear in both cases, strong NCI arise when the two counterions are close to each other. This means that a not negligible force will act in both cases just in the portion of space which lies in the middle of the two considered ionic configurations. This neighbouring effect, as confirmed also by some standard electron density calculations (not shown here), perturbs the respective molecular orbitals. This interaction leads inter alia to a weakening of the CO bond of the CH 3 OH 2 + cation which will result in the release of a water molecule and of a methenium in both the reaction paths. Of course, as the reaction 3.1 proceeds, the neutralization between the transiently formed hydride H -and the methenium CH 3 + (Fig. 3.7-b/c) will occur via a well-known electron-transfer process [START_REF] Klopman | Chemical Reactivity and the Concept of Charge-and Frontier-controlled Reactions[END_REF].

For instance, the same process takes place in the reaction 3.3 between the transient CH 3 + and the methoxide (Fig. 3.7-e/f).

The kind of analysis performed via the NCI approach cannot be easily managed by using the Fukui functions (FF), especially for the reactants (i.e., the two ionic configurations) in question.

Indeed, it is clear that in order to detect the nucleophilic (electrophilic) attack regions, it is mandatory to consider the same configuration with N + 1 (N -1) electrons. In such a case of an ionic pair configuration, it is straightforward to understand that the extra electron will accomodate into the cation, whereas when extracting an electron from the system, it will come from the anion: a situation that distorts any kind of information achievable from this analysis. Moreover, it is a common practice to use the Fukui functions/indices and some of their related local quantities as "inner" reactivity or selectivity descriptors [START_REF] Padmanabhan | Effect of Solvation on the Condensed Fukui Function and the Generalized Philicity Index[END_REF][START_REF] Chattaraj | Philicity: a Unified Treatment of Chemical Reactivity and Selectivity[END_REF][START_REF] Fuentalba | On the Condensed Fukui Function[END_REF][START_REF] Parr | Electrophilicity Index[END_REF] that are able to elucidate which specific atom will participate in a specific kind of reaction within a given molecule. To our purposes, in addition, Fukui indices are essential in the fundamental understanding of the re-distribution of the reactivity atomic centers due to the transformation of a given speciesi.e., the neutral methanol -into its cationic or anionic counterpartsi.e., CH 3 OH 2 + or CH 3 O -. In fact, the behavior of the specific relative atomic contributions to the Although it could be important to identify and visualize the spatial regions in which a specific chemical attack may occur, it is worth to be more quantitative. All the plotted spatial FF are indeed susceptible to a more or less arbitrary choise of an isosurface value. To this purpose a detailed analysis of the Fukui indices has been carried out by exploiting four different XC functionals. The first (i.e., PBE [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF]) is the XC functional employed in the Car-Parrinello Molecular Dynamics simulation described in the previous Sections. The BLYP [START_REF] Becke | Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior[END_REF] and the PW91 [START_REF] Wang | Correlation Hole of the Spin-polarized Electron Gas, with Exact Small-wave-vector and High-density Scaling[END_REF] have been tested because their respective hybrid versions are the best for this kind of calculations [START_REF] De Proft | Calculation of Molecular Electrostatic Potentials and Fukui Functions using Density Functional Methods[END_REF] while the last (i.e., PZ [START_REF] Perdew | Self-interaction Correction to Density-functional Approximations for MAny-electron Systems[END_REF]) has been used because it is builded via the Local Density Approximation; all the performed calculations show that the achievable physical insights are undependent on the nature of the XC functional used.

By evaluating the wavefunctions projections over the atomic orbitals, a Löwdin population analysis has been performed in order to obtain the condensed-to-atom Fukui indices and clarify the relative atomic contributions to the global reactivity of each of the abovementioned species.

Methanol molecule displays two relevant reaction sites for nucleophilic f + and electrophilic f -attacks, as shown in Fig. 3.9-a/b and quantitatively in Table 3.3. As expected indeed the hydrogen atom is responsible for the regional propensity to undergo an attack from a nucleophile, whereas the oxygen atom is the most susceptible site for an electrophilic attack.

This feature leads in the condensed phase to the formation of the H-bond network since the alcoholic hydrogen atoms can be regarded as the (electronic) acceptor sites and the oxygen ones as the donors.

CH 3 OH PBE BLYP PW91 PZ (LDA) Atom f + f - f 0 f + f - f 0 f + f - f 0 f + f - f 0 H 1 0.
The situation becomes different when the methanol cation is considered. As shown in Table 3.4, the positive extra-charge leads to a delocalization of the reaction centers for several type of chemical processes. Indeed the carbon atom, while almost totally inert in the neutral methanol structure, is here comparably reactive with other atomic sites. As an example, its values of the Fukui indices for nucleophilic and electrophilic attacks in principle suggest the capability of this atomic species to experience both kind of chemical reactions, especially when the PBE XC functional is employed. Moreover, it is important to remark that in this kind of calculation the neighbouring effects are not taken into account and when a nucleophile species, such as the oxygen of the methoxide, is closer to the cationic methyl group, the "nucleophilic appeal" of the carbon atom increases through polarization effects. The nucleophilic power of carbon atom relative to that of the alcoholic hydrogen(s)i.e., the f +,C /f B, the relative local philicities patterns of various atomic sites of the neutral methanol and of its cation, evaluated by employing the PBE XC functional, are as follows:

f + f - f 0 f + f - f 0 f + f - f 0 f + f - f 0 H 1 0.
• CH 3 OH -ω + : H > O > C; ω -: O > H > C; ω 0 : O > H > C. • CH 3 OH 2 + -ω + : H > C > O; ω -: O > C > H; ω 0 : O > C > H ,
where H is an alcoholic hydrogen atom. This enhancement of the carbon reactivity coupled with the homogeneization of the local chemical selectivity, when assisted by a specific local environment, can lead to the release of the CH 3 + . This latter, as a transient species, can easily combine with an H -or with a CH 3 O -, as shown in Fig. 3.7-b and Fig. 3.7-e, respectively. Of course, the enhancement of the atomic philicities is further achievable when an external electric field is applied [START_REF] Chattaraj | Chemical Reactivity Theory -A Density Functional View[END_REF].

The physical insights stemming from the investigation of the methoxide Fukui indices are informative in the same way of the FF plotted in Fig. 3.11: the oxygen atom is the most reactive site. In particular, it is of course the preferred site for an electrophilic attack. As shown also in Table 3.5, the second most reactive species, the methyl hydrogen atoms, does not show an appreciable susceptivity to nucleophilic or electrophilic attacks. Hence, here uniquely solvent and external effects are responsible for the release of the H -species when an advantageous topological arrangement of its molecular partner exists (see Fig. 3.7-a/b). The carbon atom appears to be strongly inert in this ion and the Fukui analysis fails if a LDA is employed, leading to the negative index for the nucleophilic attack f + (see the PZ (LDA) data column of Table 3.5).

Direct electric field effects can be important in evaluating the induced enhancement of the nucleophilic character at a molecular level of nucleophile species such as CH 3 O -. Indeed the responses of liquid methanol to strong electric fields are multifaceted, involving not only polarization effects and ions formation, but also a net decrease of the mean lifetime of the Hbonds [START_REF] Cassone | Liquid Methanol under a Static Electric Field[END_REF]. It is well-known that these latter in protic polar solvents stabilize the nucleophiles reducing their reactivity. Hence, methoxide, which is per se a very good nucleophile characterized by a minimum steric inhibition, gets more nucleophilic in character through a relevant increase in the H-bonds fluctuations.

Some comments within the Hard and Soft Acids and Bases (HSAB) [START_REF] Pearson | Hard and Soft Acids and Bases[END_REF] "theory" are in order.

First of all, the double channel reaction depicted in Fig. 3.7 comprises two different reaction paths which are both hard-hard in character in the first analyzed step: CH 3 OH 2 + is an hard acid and CH 3 O -is an hard base. On the other hand, the interaction between these two species, which is mostly (but not totally as we have shown through the NCI and the Fukui analysis on the CH 3 OH 2 + ) charge-controlled, splits in two reactions that are deeply different in kind.

On the one hand, it has been observed the formation of the CH 3 + and the H -as transient species, which are a soft acid and a soft base, respectively, as shown in Fig. 3.7-b. On the other hand, the formation of CH 3 + that instantaneously reacts with the methoxide (Fig. 3.7-e),

CH 3 O - PBE BLYP PW91 PZ (LDA) Atom f + f - f 0 f + f - f 0 f + f - f 0 f + f - f 0 O 0.

(Preliminary) Free Energy landscape of formaldehyde formation

As extensively treated in the previous Sections and in Ref. [START_REF] Saitta | Experiments in Atomistic Computer Simulations[END_REF], the application of a static electric field is capable to induce disparate chemical reactions and change the nature of the starting sample. In liquid methanol a field strength of 0.55 V/ Å is able to break some covalent bonds in the methyl groups of methoxide molecules and thus lead to the formation of formaldehyde, methane, and water molecules (see reaction 3.1), as copiously highlighted in this thesis. Although its microscopic mechanism and the associated atomistic reactivities have been already treated, an even more important aspect is currently under investigation: its energetics.

To the aim of investigating reaction 3.1, a series of Born-Oppenheimer simulations have been carried out in the gas phase in conjuction with a very recently presented path-CV MetD approach [START_REF] Pietrucci | Formamide Reaction Network in Gas Phase and Solution via a Unified Theoretical Approach: Toward a Reconciliation of Different Prebiotic Scenarios[END_REF]. As briefly mentioned in the Introduction and in the relative Theoretical Background Section, although this method represents an extremely powerful tool in sampling the CV space and thus in accurately reproducing the free energy landscape of very simple chemical reactions, when a relatively complex system (such as methanol) is tackled, MetD simulations based on the path-CV show problems in two-dimensional landscapes containing several (up to 5) free energy minima. This means that if one is interested only in the location of the relevant free energy basins of the reaction, the exclusive employment of MetD is sufficient whereas, when an accurate estimation of the free energy is desired, other techniques such as US [START_REF] Torrie | Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling[END_REF] can be suitably exploited in a complementary fashion.

For the gas phase variant of reaction 3.1 we limited our investigation to the exploration of the accessible basins by exploiting the just mentioned extension of the path-CV MetD [START_REF] Pietrucci | Formamide Reaction Network in Gas Phase and Solution via a Unified Theoretical Approach: Toward a Reconciliation of Different Prebiotic Scenarios[END_REF].

This preliminary study has been performed in the zero-field regime and in the presence of a moderate field (0.25 V/ Å). By comparing Figs. 3.12 and 3.13, a field-induced increase of the reaction selectivity can be straightforwardly appreciated. Indeed, while in the zero-field case the reactants (two methanol molecules) have to overcome three basins, in which the associated (meta)stable configurations are described in the caption of Fig. 3.12, in the presence of the field (Fig. 3.13) only one intermediate state survives.

Although the field-induced selective character of the gas phase reaction is an interesting indication of the role played by an electric field in assisting a given chemical reaction, more exhaustive and quantitative investigations are currently under study in the condensed phase at several field strengths. Indeed the solvent appears to be essential in lowering the relative energy barrier between the initial and the final states. Also in the liquid state, an "exploratory"

MetD calculation in the zero-field regime with the extended path-CV MetD approach has been already performed. It revealed that only the reactants and the products basins are relevant (i.e., the reaction occurs directly without involving other intermediate metastable states). Also in the presence of the electric field we observed that it was possible to converge the underlying probability distributions in a relatively confined region of the CV space, in which only the products and the reactants states were present.

In order to achieve statistical convergence, the final sampling on the CV space has been conducted via the US [START_REF] Torrie | Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling[END_REF] technique, as previously mentioned. In particular, only the zero-field regime "umbrellas" have (at this moment) a reliable statistical convergence. For this particular case, 66 parabolic potentials have been employed in order to fill and probe the CV space and, in some critical cases, also in conjuction with proper linear corrections to the functional form of the bias potential. These 66 independent biased simulations were recollected by means of WHAM [START_REF] Kumar | The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method[END_REF][START_REF] Ferrenberg | Optimized Monte Carlo Data Analysis[END_REF][START_REF] Bennet | Efficient Estimation of Free Energy Differences from Monte Carlo Data[END_REF][START_REF] Shirts | Statistically Optimal Analysis of Samples from Multiple Equilibrium States[END_REF], which is capable to reconstruct the underlying free energy, as shown in Fig. 3.14. As it is clear, the energy barrier that the reactants have to overcome without the electrical assistance is of order of 50 kcal•mol -1 , making impracticable the evolution of the reaction.

The preliminary results on the free energy landscape of this chemical reaction in the presence of an electric field (at various strengths) suggest that this latter is capable not only to drive the reaction via a strong reduction of the energy barrier, but also to "open" and/or "close" peculiar reaction pathways (i.e., the transition states) as succinctly shown in Table 3.6. The statistically converged evaluation of the FESs at different field strengths is currently under investigation. Finally, the solvent plays a central role in assisting this chemical reaction as it has been possible to observe via an indicative comparison of the free energies associated with the gas phase and the liquid phase cases. Further analysis are ongoing.

circumstance that prevents the transfer of a proton along the same path.

The investigation of another sample led us to highlight some relevant properties in assisting the molecular dissociation and the ionic conduction. In this work we thus analyzed the effects of a static and spatially uniform electric field on bulk liquid methanol.

We found that the effect of the field on the local structure of the liquid, as monitored through the radial correlation functions, critically depends on its intensity: in fact, weak and moderately intense fields stabilize the H-bond network while too strong fields eventually break the bonds and disrupt the network. Correspondingly, the molecular dipole moment gets, on average, more and more aligned with the field direction, a circumstance that favors the transfer of protons along the H-bond structure until molecular dissociation occurs for a field strength of about 0.31 V/ Å, a value slightly larger than that for water. However, we observed a net proton flow for fields above 0.36 V/ Å with an ohmic current-voltage diagram as is also the case for water.

Although it is evident that the underlying mechanism of ionic defects migration is Grotthusslike, the proton transport efficiency is at least one order of magnitude lower in methanol as compared with its counterpart for water. We interpreted this difference by considering the field-induced fluctuations in the H-bond distribution of the ionic solvation shells. In fact, these fluctuations enhance the matching probability between the coordination number of a charge carrier such as an hydronium ion and that of a neutral water molecule, thus making the ionic transfer process possible. However, in liquid methanol the coordination numbers of a neutral CH 3 OH molecule and that of the ionic species CH 3 OH 2 + are the same, a circumstance that does not produce an enhancement of the proton transfer efficiency via H-bond fluctuations as in liquid water.

Finally, above the field threshold of about 0.55 V/ Å, the chemical nature of a sample of pure liquid methanol dramatically changes. Indeed, the formation of several chemical species has been observed starting from this field intensity. Among the others, the most important formed compounds are formaldehyde (H 2 CO), methane (CH 4 ), and dimethyl ether (CH 3 OCH 3 ). Con-cerning the previously unknown chemical reactions in which the simultaneous formation of H 2 CO and CH 4 has been observed, a preliminary investigation of the associated energetics, via the evaluation of the free energy surface in the zero-field regime and only partially also in presence of electrical perturbations, has been performed. This preliminary study revealed that the electric field strongly affects the evolution of the reaction by means of not only a drastic change of the mutual stability of the reactants and the products basins but also of the accessible reaction pathways.

By comparing this reaction mechanism with that leading to the formation of CH 3 OCH 3 , a double channel reaction can be argued. In particular, in both reactions the putative transition state is the ionic pair composed by methoxide CH 3 O -and the methanol cation CH 3 OH 2 + .

Moreover, depending on their mutual orientation and independently on their orientation with respect to the field direction, a chemical reaction can proceed leading to H 2 CO and CH 4 , in one case, and to CH 3 OCH 3 in the other. In order to disclose the relative atomistic reactivities, a detailed study of the Fukui indices has been thus performed. This latter revealed an increased general reactivity of the carbon atom when considered within the cationic state CH 3 OH 2 + rather than in the neutral methanol CH 3 OH one. This finding, among the others, accounts for the observed chemistry of methanol under the effect of a static and homogeneous electric field.

Since chemical reactions occur almost always between different molecular species it is important to know in which proportions the two molecules are somewhat responsible for what is happening. Recently, from an empirical relationship [START_REF] Maynard | Reactivity of the HIV-1 Nucleocapsid Protein p7 Zinc Finger Domains from the Perspective of Density-functional Theory[END_REF] an electrophilicity index ω has been developed [START_REF] Parr | Electrophilicity Index[END_REF]. This parameter can be related to the others because it is possible to expand where the chemical hardness η appears. The role played by this latter (or by the chemical softness S, obviously) is not negligible as it enters as a protagonist in the HSAB principle [START_REF] Parr | Absolute Hardness: Companion Parameter to Absolute Electronegativity[END_REF].

Therefore also the function f (r)i.e., the Fukui function -must hold a relevant part and by exploiting the Maxwell relation it is clear that this function is defined as the differential change in electron density due to an infinitesimal change in the number of electrons:

f (r) =   ∂n(r) ∂N   v(r) . (A.11)
Because of the derivative discontinuity for isolated molecules at zero temperature, the FF are defined from above and below using the one-sided derivatives and in such a way that one has f + (r) and f -(r), respectively. For practical purposes and since the number of electrons can vary, in principle, in a discrete manner, the FF can be simply written as .13) which in the frozen orbital approximation can be seen as the squares of LUMO and HOMO (i.e., their density). Because of the FMO regioselectivity analogs that can be made for nucleophilic (f + v,N ) and electrophilic (f - v,N ) attack, it is possible to define a reactivity indicator for a radical attack:

f + v,N (r) = n v,N +1 (r) -n v,N (r) (A.12) f - v,N (r) = n v,N (r) -n v,N -1 (r) (A
f 0 v,N (r) = 1 2 (f + v,N (r) + f - v,N (r)) = 1 2 (n v,N +1 (r) -n v,N -1 (r)). (A.14)
From the FF denition stems the fundamental concept for which they can be regarded as an important quantity in identifying the chemical character of a given portion of space in a molecule. Indeed it can be demonstrated that for a given amount of charge transfer ∆N , the most stable way to distribute the additional charge in a molecular species corresponds to ∆n(r) = ∆N f (r) [START_REF] Ayers | Variational Principles for Describing Chemical Reactions: the Fukui Function and Chemical Hardness Revisited[END_REF]. As shown, although the name of these functions has been taken from the pioneer of the frontier orbital reactivity descriptors, the concept of FF is entirely DFTbased; therefore in principle the achievable results with them are correct, taking into account the electron correlations and orbitals relaxations.

A.2 Local reactivity descriptors

Although the global FF have their importance, it is more likely that one is interested not in having the reactivity description of a point in space but rather to give a reactivity indicator to each atom in a molecule. More precisely, one whishes to identify which atom in a molecule is most likely to react with an attacking electrophile or nucleophile. To this purpose it has been developed the concept of condensed-to-atom Fukui in which it is possible to have an expression depending on the atomic populations [START_REF] Yang | The Use of Global and Local Molecular Parameters for the Analysis of the Gas-Phase Basicity of Amines[END_REF]: or, in terms of the atomic charges q (α) ,

f +,α v,
f +,α v,N = q (α) v,N +1 -q (α) v,N
(A.17)

f -,α v,N = q (α) v,N -q (α) v,N -1 (A.18) f 0,α v,N = 1 2 (q (α) v,N +1 -q (α)
v,N -1 ) , (A. [START_REF] Tajima | Phase Transition in KOH-Doped Hexagonal Ice[END_REF] where α labels the atom and these relationships are for the nucleophilic, electrophilic, and radical attack, respectively.

The condensed-to-atom Fukui concept is extremely important in chemical DFT because it represents the link between the global and the local reactivity descriptors. Indeed it is possible to reformulate the theoretical framework presented in the previous Section with an atomic resolution by modulating each global reactivity indicator via the condensed Fukui (i.e., Eqs. A.17, A.18, A.19). After having chosen a given population analysis scheme such as that by Hirshfeld [START_REF] Hirshfeld | Bonded-Atom Fragments for Describing Molecular Charge Densities[END_REF] or Löwdin, is then possible to easily evaluate the condensed-to-atom softness from the following relation

s k,α = S • f k,α (k = +, -, or 0) (A.20)
which, obeying to the HSAB principle, can be exploited in order to predict which atom of a nucleophile species will interact with which atom of an electrophile molecule [START_REF] Chattaraj | Chemical Reactivity Theory -A Density Functional View[END_REF]. Indeed it is possible to show that considering the ith atom of a species A and the jth atom of a species B, then the most favorable situation that arises from the minimization of the grand canonical potentials of A and B leads to s Ai =s Bj .

Up to now the introduced mathematical objects are of a certain importance which is unfortunately constrained by the fact that FF and their condensed-to-atom equivalents are always referred to the given molecule to which they refer. Indeed they are defined at fixed external potential and no comparisons can be made between different molecular entities. On the other hand, it has been recently introduced a local quantity called philicity [START_REF] Chattaraj | Philicity: a Unified Treatment of Chemical Reactivity and Selectivity[END_REF] which is the local counterpart of the electrophilic power ω defined in Eq. A.7 and which need again the condensed-to-atom FF aid:

ω k,α = ω • f k,α (k = +, -, or 0) . (A.21)
This local reactivity descriptor, as well as the local softness, have the peculiarity to take advantage from the local regiosensitivity which is proper of the condensed-to-atom FF and at the same time to allow to perform reactivity comparisons between a set of (not so) different molecules. In particular, it has been recently shown that in a series of substitued molecules, when the substituent represents a small portion of the molecule and the softnesses are very similar, the philicity can give a better resolution of the sites reactivity of the species [START_REF] Cedillo | A Local Extension of the Electrophilicity Index Concept[END_REF]. The last selectivity index having an atomic resolution is the local hardness [START_REF] Meneses | A Proposal for a New Local Hardness as Selectivity Index[END_REF]. This condensed-toatom reactivity descriptor is again related to the definition of FF but encompasses the necessity of imposing the knowledge of having an electrophilic, a nucleophilic, or a radical attack. Making use of the HOMO and LUMO energies of the neutral molecule, this descriptor is proper to describe the intramolecular selectivity regions and it is defined as in the following:

η α = ǫ LU M O f +,α -ǫ HOM O f -,α , (A.22)
where ǫ LU M O and ǫ HOM O are the LUMO and HOMO energies, respectively. The local hardness reduced density gradient, stemming from the electron density itself and from its first derivative, s = 1 2(3π 2 ) 1/3 |∇n| n 4/3 , (A. [START_REF] Cwiklik | Hydroxide Trapped in the Interior of Ice: a Computational Study[END_REF] is also a key quantity since it is employed in order to describe the deviation from a homogeneous electron distribution [START_REF] Hohenberg | Inhomogeneous electron gas[END_REF][START_REF] Becke | Modern Electronic Structure Theory[END_REF][START_REF] Cohen | Insights into Current Limitations of Density Functional Theory[END_REF]. Moreover, its properties were a crucial aspect during the development of new and more accurate functionals [START_REF] Zupan | Distributions and Averages of Electron Density Parameters: Explining the Effects of Gradient Corrections[END_REF].

The NCI (or noncovalent contacts) can be identified with the regions of the real space characterized by a small reduced density gradient at low electron density [START_REF] Contreras-García | NCIPLOT: a Program for Plotting Noncovalent Interaction Regions[END_REF][START_REF] Johnson | Revealing Noncovalent Interactions[END_REF]. In their implementation for the visualization, these regions are mapped by plotting an isosurface of s for a reasonably low value of the reduced density gradient. Moreover, the sign of the second eigenvalue of the Hessian times the density itself can be used in order to identify different regions in the isosurfaces maps allowing to discerne the intensity as well as the stability of the nature of these interactions. At the same time, the sign of the Laplacian ∇ 2 n gives informations about the net flux, i.e., ∇ 2 n < 0 (near the nuclei) if it is incoming and ∇ 2 n > 0 if it is outcoming from an infinitesimal volume in the neighbourhood of a given reference point.

In order to distinguish each weak interaction, the contributions to the Laplacian along the axes of its maximum variation have to be analyzed. These are the eigenvalues λ i of the Hessian matrix of the electron density (second derivatives) such that

∇ 2 n = λ 1 + λ 2 + λ 3 , (λ 1 < λ 2 < λ 3 ) . (A.24)
In the neighbourhood of the nuclei, all the λ i are negative, whereas far from them λ 3 > 0.

Obviously, in molecules the λ 3 values can vary along the internuclear direction while λ 1 and λ 2 gives informations about the fluctuations of n in the orthogonal plane of the eigenvector of λ 3 ; λ 2 can be either positive or negative depending on the nature of the interaction. Bonding interactions such as the H-bond are characterized by an accumulation of electron density

Introduction 1 1 . 2

 112 Theoretical Background 1.1 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Hohenberg-Kohn Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 The Levy approach and the Kohn-Sham scheme . . . . . . . . . . . . . . 1.1.3 The exchange and correlation functional . . . . . . . . . . . . . . . . . . 1.1.4 Plane waves and pseudopotentials . . . . . . . . . . . . . . . . . . . . . . Ab initio simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Car-Parrinello Molecular Dynamics . . . . . . . . . . . . . . . . . . . . . 1.2.2 Static electric fields in ab initio simulations . . . . . . . . . . . . . . . . . 1.3 Collective variables and Metadynamics . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Path Collective Variables with a new definition of distance . . . . . . . . 1.3.2 Umbrella Sampling and the Weighted (and Dynamic) Histogram Analysis Methods . . . . . . . . . . . . . . . . 2 Results on Proton Transfer in H-bonded Systems 2.1 Water ices under an electric field . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Liquid methanol under an electric field . . . . . . . . . . . . . . . . . . . . . . . 3 Results on Methanol Chemical Reactions 3.1 Methanol chemistry at high electric fields . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Field-induced chemical reactions . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Reactivity and (regio)selectivity . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 (Preliminary) Free Energy landscape of formaldehyde formation . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusions A Theoretical background on Fukui Functions, (Regio)Selectivity, and Reactivity Descriptors A.1 Global reactivity descriptors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Local reactivity descriptors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Another approach: the Noncovalent Interactions . . . . . . . . . . . . . . . . . . B Supplemental Results on the Reactivity Descriptors B.1 Global and (other) local indicators . . . . . . . . . . . . . . . . . . . . . . . . .

(1. 45 )

 45 t w is the kinetic energy density, and b 1 = 0.04918, b 2 = 0.132, b 3 = 0.2533, and b 4 = 0.349.

Figure 1 . 1 :

 11 Figure 1.1: Construction principle of topological path-CV [36]. Arrows highlight changes of coordination numbers whereas all other matrix elements are free to change. Courtesy of the authors of Ref. [36].

Figure 2 . 1 :

 21 Figure 2.1: Relative energies corresponding to proton shifts of a given distance along an H-bonded chain formed by three water molecules, for an electric field intensity of 0.36 V/ Å, in a system with 64 molecules; the solid blue line refers to ice XI while the other curves refer to ice I h and were generated using three different molecular configurations that are shown in the inset; panel (a): configuration of three water molecules present in both I h and XI structures leading to the blue curve and green curve with circles, respectively, upon displacing protons as discussed in the text; panel (b): configuration of three water molecules in ice I h leading to the red line with crosses; panel (c): configuration of three water molecules in ice I h leading to the black line with diamonds.

Figure 2 . 2 :

 22 Figure 2.2: Donor-acceptor distances during the first ionization events (identified by the arrows) in a system with 64 molecules: the solid blue line refers to ice XI while the dashed red line refers to ice I h . The inset shows (a) the standard ice XI structure and (b) the field-induced planes flattening, with pertinent angles labeled as Greek letters; panel (a): α = β = γ = δ = ǫ = ζ = 109.45 • ; panel (b): α = 113.2 • , β = 112.9 • , γ = 102.3 • , δ = 114.8 • , ǫ = 114.8 • , and ζ = 102.4 • , a configuration which exhibits a shortening of the in-plane intermolecular lengths.

Figure 2 . 3 :

 23 Figure 2.3: Ionic current-voltage characteristic for an orthorhombic cell with edges a = 9.05 Å, b = 15.67 Å and c = 14.87 Å: solid circles: samples of ice I h with 64 molecules; open circles: samples of ice XI with 64 molecules; diamonds: samples of ice I h with 96 molecules; triangle: samples of ice I h with 128 molecules (lines are a guide for the eye).

  order to evaluate the equation of state (EOS), we perform a series of relaxation simulations by applying the field along the c-axis in ices I h and XI, and along the b-axis in the protonordered structure only. In this way we obtain the energy profiles versus deformation over the range of the explored field intensities that are shown in Fig. 2.4. In ice XI, when the field is oriented along the c-axis, upon increasing the field strength the minimum of the energy shifts to lower values of the ratio between the cell parameters c and a (Fig. 3.4-a); instead, no shift at all is observed in the I h phase (Fig. 3.4-b) or when the field is oriented along the b-axis of the ferroelectric sample (Fig. 3.4-

Figure 2 .Figure 2 . 4 :Figure 2 . 5 :

 22425 Figure2.5 shows the elastic constant c 66 , plotted as a function of the field intensity; in ice I h the value of this quantity computed for vanishing electric field (1.9 GPa) turns out to be lower than both the available experimental results[START_REF] Petrenko | Physics of Ice[END_REF][START_REF] Fletcher | The chemical physics of ice[END_REF][START_REF] Gammon | Elastic Constants of Artificial and Natural Ice Samples by Brillouin Spectroscopy[END_REF] and another independent theoretical estimate[START_REF] Penny | A Theoretical Determination of the Elastic Constants of Ice[END_REF] by about 1.5 GPa. A manifestly different behavior is observed in the two crystalline structures: in ice I h the elastic constant c 66 shows a weak dependence on the field strength, while in ice XI it actually decreases with the field intensity in a linear way. This difference is ultimately related to the structural symmetries. The proton-ordered sublattice of ice XI allows for an enhanced coupling between the molecular dipole moments and the external field when this is oriented along the ferroelectric axis. It turns out that this circumstance makes the system softer by forcing a partial molecular alignment along the field direction. In the cases of ice I h subject to a field along the c-axis and of ice XI under the effect of an electric field acting along the b-cell parameter direction, the coupling of the local molecular dipole moments with the external field is not as straightforward; the different distributions of dipole moments,

Figure 2 . 6 :

 26 Figure 2.6: Donor-acceptor distance between two molecules measured on the occasion of the first molecular dissociation events when the field is oriented along the b-axis of ice I h (red dashed line) and of ice XI (blue solid line). The instant in which the ionization occurs is signalled by the arrow.

Figure 2 . 7 :

 27 Figure 2.7: a) Average donor-acceptor distances of an H-bond chain of ice I h in which the first dissociation

Figure 2 .

 2 Figure 2.8-b/c shows the comparison with the previously investigated cases where the fieldis oriented along the c-axis for ice I h and XI, respectively[START_REF] Cassone | Proton Conduction in Water Ices under an Electric Field[END_REF]. All the ionic current-voltage diagrams show a ohmic behavior after a net PT process has taken place. These results confirm those found in several experiments performed on ice I h samples[START_REF] Maidique | Transfer of protons through pure ice in single crystals -Part III. The dielectric relaxation spectra of water, ice, and aqueous solutions and their interpretation[END_REF][START_REF] Liu | Determination of Phase Diagram of Water and Investigation of Electrical Transport Properties of Ices VI and VII[END_REF]. It is not surprising to observe a good agreement between the two different situations (see Fig.2.8-b); a similar conclusion follows from the calculations performed with the BLYP functional. However, as far as ice XI is concerned, we observe a different behavior below the current threshold value (Fig.2.8-c). This result can be interpreted at relatively low voltages (up to about 5 V for

Figure 2 . 8 :

 28 Figure 2.8: a) Ionic current-voltage diagrams corresponding to a field oriented along the cell parameter b=15.67 Å; red solid diamonds refer to a sample of ice I h ; blue solid circles refer to a sample of ice XI. b) Ionic current-voltage diagrams of ice I h under the action of a field oriented along the b-axis (red open circles) and the c-axis (red solid circles). c) Ionic current-voltage diagrams of ice XI under the action of a field oriented along the b-axis (blue open circles) and the c-axis (blue solid circles).

Figure 2 . 9 :

 29 Figure 2.9: The Cmc2 1 ice XI orthorhombic centered unit cell (partially modified from Fig. 1 in Ref. [20]).

( 3 -

 3 and 4-type in Fig.2.9) are hosted which exhibit a large positive y-component of dipole moments. In the absence of the field, this component falls in a range of 2.51 D-2.69 D and 2.90 D-3.10 D for the 3-and 4-type molecules shown in Fig.2.9, respectively[START_REF] Batista | Molecular Multipole Moments of Water Molecules in Ice I h[END_REF][START_REF] Ishii | First-Principles Study of Spontaneous Polarization and Water Dipole Moment in Ferroelectric Ice XI[END_REF]. One half of the ice XI planes are ferroelectric with respect to the b-axis direction and the remaining ones are anti-ferroelectric.

  2.10). Moreover, the condition of partial ferroelectricity is only a non-equilibrium situation induced by the field and by the related conspicuous amount of defects (seeFig. 2.11). Because of the formation of defects and their propagation via field interactions, one can argue that the conduction properties of ice XI are very similar to those discussed in the two cases in which a field along the ferroelectric axis or along the b one is involved, provided that the current threshold of the latter has been reached and overcome.

Figure 2 . 10 :

 210 Figure 2.10: "Impossible" ferroelectric phase (along the b-axis) of ice XI.

Figure 2 . 11 :

 211 Figure 2.11: A realistic sketch, without H-bond perturbations, of the molecular arrangement when an electric field stronger than 0.29 V/ Å is applied along the b-axis of ice XI. Some molecules are not shown.

  way H-bonds model the structure of liquid methanol is manifest in the oxygen-oxygen (O-O) and oxygen-hydrogen (O-H) radial distribution functions (RDF) that are shown in Fig. 2.12; carbon-carbon (C-C), carbon-oxygen (C-O), carbon-hydrogen (C-H), and hydrogen-hydrogen (H-H) RDFs are also plotted in Fig. 2.13. As usual, hydrogen atoms involved in the computed RDFs are those bonded to the oxygen atom in the methanol molecule.

Figure 2 . 12 :

 212 Figure 2.12: Oxygen-oxygen (a) and oxygen-hydrogen (b) radial distribution functions (continuous lines) plotted, with pertaining coordination numbers (dashed lines), as a function of the radial separation for zero electric field.

Figure 2 . 13 :

 213 Figure 2.13: Carbon-carbon (a), carbon-oxygen (b), carbon-hydrogen (c), and hydrogen-hydrogen (d) radial distribution functions plotted as a function of the radial separation for zero electric field.

Figure 2 . 14 :

 214 Figure 2.14: Oxygen-oxygen radial distribution function calculated for increasing values of the electric field. Black solid line: E = 0 V/ Å; red dashed line: E = 0.15 V/ Å; green dotted line: E = 0.31 V/ Å; blue dash-dotted line: E = 0.46 V/ Å. Inset: heights of the first peak (black solid circles) and of the first dip (red solid circles) plotted as a function of the field intensity.

Figure 2 . 15 :

 215 Figure 2.15: Oxygen-hydrogen radial distribution function calculated for increasing values of the electric field. Black solid line: E = 0 V/ Å; red dashed line: E = 0.15 V/ Å; green dotted line: E = 0.31 V/ Å; blue dash-dotted line: E = 0.46 V/ Å. Inset: heights of the first peak (black solid circles) and of the first dip (red solid circles) plotted as a function of the field intensity.

Figure 2 . 16 :

 216 Figure 2.16: O• • •HO angle distribution function plotted for several electric field strengths. Black solid line: E = 0 V/ Å; red dashed line: E = 0.10 V/ Å; green dotted line: E = 0.26 V/ Å; blue dash-dotted line: E = 0.36 V/ Å. The peak of the distribution function shifts to lower values of the O• • •HO angle as the field strength increases up to 0.26 V/ Å (green dotted line).

Figure 2 .

 2 Figure 2.16 shows the ∠HO• • • O bond angle distribution function: for increasing values of the field intensity the distribution gets sharper, the height of the peak increases while the position shifts to lower values of the angle. A similar trend was found by Sun et al. [153] in a classicalMolecular Dynamics which, however, did not capture the reentrant feature discovered with the ab initio scheme. As noted before for other structural properties of methanol, the reported trend is reversed for fields approximately stronger than 0.3 V/ Å.As far as the other structural properties are concerned, we found that the C-C and C-O RDFs of liquid methanol do not show a significant sensitivity to the external field intensity, at least for not too strong fields; just the C-H RDF exhibits a weak dependence on the field strength similar to that observed in the O-H RDF.

Figure 2 . 17 :

 217 Figure 2.17: Distribution of the angle θ formed by the molecular dipole moment vector with the electric field (z) axis. Black curve: E = 0; red curve: E = 0.26 V/ Å.

Figure 2 . 18 :

 218 Figure 2.18: Time autocorrelation function of H-bonds for the zero-field case (black solid line) and for a field strength of 0.36 V/ Å (red dotted curve).

Figure 2 .

 2 Figure 2.19: "Zundel-like" ionic complex [C 2 H 6 O 2 H 3 ] + during a PT event.

Fig. 2 .Figure 2 . 20 :

 2220 Fig. 2.19) forms during the PT process, but the molecular cooperativity turns out to be much lower than in liquid water, presumably because of steric effects; the [C 2 H 6 O 2 H 3 ] + complex can be considered as the counterpart of the Zundel ion that is involved in the PT mechanism in liquid water. As a result of the molecular dissociation induced by the electric field, the sample is characterized by the presence of ionic species such as CH 3 OH 2 + and CH 3 O -. The local H-bond environment is obviously perturbed by charge defects; at the same time, as already discussed, the H-bond lifetime decreases. These two circumstances lead in liquid water to an enhanced PT probability as a result of the increased fluctuations of the first solvation shell of H 3 O + species. In this way, the hydronium ion can propagate, with corresponding variations of the coordination numbers of the solvation shells. In fact, the first coordination shell of H 3 O + includes, on average, three

Figure 3 . 1 :

 31 Figure 3.1: Carbon-carbon radial distribution function determined at different field strengths. At and above a threshold of 0.55 V/ Å the system shows the onset of new chemical species. Pink and light green triangles are referred to field intensities of 0.55 V/ Å and 0.60 V/ Å, respectively.

  parts. If in the system there is a sequence of adjacent molecules (and ions) such as H 2 O, H 2 CO, CH 3 OH, and CH 3 O -, as schematically shown in the left side of Fig.3.2, the water molecule can dissociate under the field action and hence donate a proton to its H-bonded formaldehyde which obviously becomes an H 2 COH + cation (Fig.3.2, central sketch). In concert, the nearest methanol can lose a proton in favor of the adjacent methoxide along the field direction. The neutralization of the two resulting ions, H 2 COH + and CH 3 O -, is the last natural consequence of the effect of the applied field that gives rise to the methoxymethanol formation (Fig.3.2, right).

Figure 3 . 2 :

 32 Figure 3.2: Methoxymethanol formation mechanism. Starting from water, formaldehyde, methanol, and the anion of this latter (left), the displacement of the negative charge against the electric field direction (center and right) leads to a recombination process which produces the methoxymethanol molecule (right).

Figure 3 . 3 :

 33 Figure 3.3: Dimethyl ether and the least favored formaldehyde formation mechanisms. a) Starting from a

Figure 3 . 4 :

 34 Figure 3.4: Dimethyl ether formation mechanism. Starting from the two common methanol ions (left), the neutralization process begins with the field-induced displacement of the positive charge toward the cation body reaching its own methyl group. The presence of the negative head of the anion not only enhances this process but also contributes to the release of a methenium (center) which can bonds with it giving rise to dimethyl ether and water (right).

Figure 3 . 5 :

 35 Figure 3.5: Formaldehyde and methane formation mechanism. Starting from the two common methanolions (left), the neutralization process begins with the field-induced displacement of the positive (negative) charge along the cation (anion) body reaching its own methyl group. The charge recombination occurs via a methenium and an hydride release stemming from the cation and the anion, respectively (center). This leads to the formation of water, methane, and formaldehyde (right).

Figure 3 . 6 :

 36 Figure 3.6: Concerted cooperative molecular motions which allows the onset of a formaldehyde in the sample (a). This species will be employed for the formation of a dimethyl ether molecule (b). These sketches show also the key role played by the solvent in assisting these chemical reactions as a reservoir of proton and hydride accepting/donating sites (b). The overall process occurs in about 300 fs.

Fig. 3 . 6 .

 36 Fig. 3.6. Here, a newly formed formaldehyde through the common process (Fig. 3.6-a) can interact with a just created methanol cation giving rise to the onset of a dimethyl ether species (Fig. 3.6-b).

3 . 3 Figure 3 . 7 :

 3337 Figure 3.7: A double channel reaction. Depending on the relative orientation of the methoxide CH 3 O -with respect to its cationic counterpart CH 3 OH 2+ , the reaction proceeds yielding formaldehyde, methane, and water (a,b,c) or dimethyl ether and water (d,e,f). These two reactions occur in the bulk liquid under the effect of an external electric field and also in the gas phase without any external perturbation with the same mechanisms (i.e., the same transition states).

Field

Figure 3 . 8 :

 38 Figure 3.8: NCIs between the two local reactants just before the formation of formaldehyde, methane, and water (a) and dimethyl ether and water (b).

Figure 3 . 9 :

 39 Figure 3.9: Fukui functions of methanol for the nucleophilic (a) and eletrophilic (b) attacks.

Figure 3 . 10 :

 310 Figure 3.10: Fukui function for the nucleophilic attack of the methanol cation.

Figure 3 . 11 :

 311 Figure 3.11: Fukui function for the electrophilic attack of the methoxide.

Figure 3 . 12 :

 312 Figure3.12: "Exploratory" (not converged: the depth of the basins is meaningless!) FES projected onto the path-CV space of reaction 3.1 performed in the gas phase and in the zero-field regime. Three basins separate the reactants and the products. In ascending order of S (progress along the reaction path) we have: 1) reactants: two methanol molecules; 2) a methanol, a formaldehyde, and a hydrogen molecule; 3) a formaldehyde monohydrate and a methane molecule; 4) a methane, a hydrogen, and a formic acid molecule; 5) products: a formaldehdye, a methane, and a water molecule.

Figure 3 . 13 :

 313 Figure 3.13: "Exploratory" (not converged: the depth of the basins is meaningless!) FES of reaction 3.1performed in the gas phase under the effect of a field strength of 0.25 V/ Å. Two of the three basins previously observed in the absence of the electrical perturbation are not accessible under these circumstances. In ascending order of S we have: 1) reactants: two methanol molecules; 2) a formaldehyde monohydrate and a methane molecule; 3) products: a formaldehdye, a methane, and a water molecule.

 9 )

 9 in a Taylor series up to the second order the energy change as follows: ∆E = E(N + ∆N, v) -E(N, v) ≈ µ∆N + 1 2 η(∆N ) 2 , (A.6) since µ and η are defined as in eqs. A.1 and A.2. If we define ∆N max as the maximum number of electrons that an electrophile can accept and ω = -∆E min as the related stabilization energy, ∆N max = -µS → ω = µ 2 2η . (A.7) Although the electrophilic power ω is an extremely important quantity, it is not able by itself to perform reactivity comparisons between different molecular sites and for this purpose a local version of this parameter has been proposed [172]. This local index distributes ω over a given molecule making use of the Fukui functions (FF). The FF are so called in honor of the Nobel Prize for Chemistry of 1981 Kenichi Fukui who was awarded for his developments of Frontier Molecular Orbital (FMO) theory "[...] concerning the course of a chemical reaction". The FF have been introduced in 1984 by Parr and Yang [29] as a consequence of other works previously published by Parr et al. [30] concerning the chemical DFT and on the already briefly treated concepts such as the chemical potential µ and the hardness η. Historically the equalization of the electronic chemical potential to the electronegativity changed in sign suggested the possibility of treating electronic structure theory in a thermodynamics fashion. Ergo, Parr et al. wrote the total differential of the energy E[N ,v(δv(r)dr = µ v,N dN + n v,N (r)δv(r)dr. (A.8) By extending the thermodynamic analogy it is possible to obtain the DFT-equivalent of the Legendre-transform structure of classical thermodynamics and the following Maxwell rela-In analogy to Eq. A.8 the total differential for the chemical potential µ[N, v(r)] is dµ = )dr = ηdN + f (r)δv(r)dr (A.10)

Figure A. 1 :

 1 Figure A.1: A typical methanol-methanol H-bond as seen in the Noncovalent Interaction formalism.

  ′ ext = v ext + cost) and the set G of resulting ground states such that |Ψ ′ 0 = e iϕ |Ψ 0 , where ϕ is any global phase. Let define the map from V to G as A:

The many-body eigenstates |Ψ k corresponding to the Hamiltonian (1.1) are obtained as a solution of the stationary Schrödinger equation; more precisely, this is true for the solution referred to the ground state |Ψ 0 on which we will focus. Let consider now the set of all Hamiltonians of the form (1.1) with non-degenerate ground state, such as the set of all the local potentials v ext for which the Schrödinger equation leads to non-degenerate ground states, whereas the interaction (1.2) is kept fixed. This set contains an infinite number of fictitious copies (i.e., an ensemble) of any given v ext (r), which are obtained by simply adding an arbitrarily constant to v ext (r). This treatment can be otherwise generalized to the case of degenerate ground states. The solution to the Schrödinger equation can be interpreted as a relation between the set V of external potentials which differ from each other for more than a constant (i.e., v

Table 2 . 1 :

 21 Percent fractions of ionic defects induced in ice I h and ice XI by increasing electric-field strengths.

	I h	Ice XI

Table 2 . 2 :

 22 Percent fractions of defects present in ice I h doped with increasing molarities of CaCl 2[START_REF] Grimm | Low-Frequency Electrical Properties of Polycrystalline Saline Ice and Salt Hydrates[END_REF].

Table 2 . 3 :

 23 Positions ( Å) of the first peak of atomic radial distribution functions for zero electric field. Column

	R 1st peak	PBE	Experiment	BLYP	MC	MD
	O-O	2.71	2.74	2.94	2.88	2.81
	C-O	3.54	3.58	3.76	3.62	3.56
	C-C	4.23	3.97	4.22	4.10	4.21
	O-H	1.71	1.77	1.91	1.96	1.87
	C-H	2.66	2.68	2.87	2.90	N.A.
	H-H	2.34	2.42	2.51	2.50	N.A.

2: present estimates obtained with the PBE exchange and correlation functional; column 3: experimental data

[START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF][START_REF] Yamaguchi | The Structure of Liquid Methanol Revisited: a Neutron Diffraction Experiment at -80 • C and +25[END_REF]

; column 4: results obtained with the BLYP functional

[START_REF] Morrone | Ab initio Molecular Dynamics Study of Proton Mobility in Liquid Methanol[END_REF]

; column 5: results obtained with Monte Carlo simulations

[START_REF] Shilov | Long Range Order and Hydrogen Bonding in Liquid Methanol: A Monte Carlo Simulation[END_REF]

; column 6: data extracted from classical Molecular Dynamics simulations

[START_REF] Sun | Molecular Dynamics Simulation of Liquid Methanol under the Influence of an External Electric Field[END_REF]

.

Table 3 . 2 :

 32 Formaldehyde mean lifetimes at different field intensities.

  2 CO, CH 4 , and H 2 O (reaction 3.1) and CH 3 OCH 3 and H 2 O (reaction 3.3), it is important to clarify the roles played by each atom. The two species that

Table 3 . 3 :

 33 Condensed-to-atom Fukui indices for nucleophilic (f + ), electrophilic (f -), and radical attack (f 0 ) for the methanol molecule evaluated with different XC functionals. In red are highlighted the highest values of the indices for each specific kind of chemical attack. Note that the atom H 1 represents the alcoholic hydrogen.

		0197 0.0916 0.0557 0.0340 0.0961 0.0651 0.0248 0.0903 0.0575 0.0430 0.0875 0.0653
	O	0.0098 0.4609 0.2354 0.0021 0.4598 0.2305 0.0117 0.4597 0.2357 0.0144 0.4579 0.2362
	C	0.0070 0.0848 0.0459 0.0177 0.0729 0.0453 0.0063 0.0758 0.0411 0.0119 0.0814 0.0933
	H 3	0.0058 0.1613 0.0836 0.0087 0.1314 0.0701 0.0085 0.1624 0.0855 0.0127 0.1626 0.0877
	H 4	0.0097 0.0882 0.0490 0.0097 0.1251 0.0674 0.0112 0.0924 0.0518 0.0141 0.0924 0.0666
	H 5	0.0058 0.1615 0.0837 0.0088 0.1594 0.0841 0.0083 0.1627 0.0855 0.0127 0.1629 0.0878

Table 3 . 4 :

 34 Condensed-to-atom Fukui indices for nucleophilic (f + ), electrophilic (f -), and radical attack (f 0 ) for the methanol cation evaluated with different XC functionals. In red are highlighted the highest values of the indices for each specific kind of chemical attack: it is visible a spread of the reactivity over the atoms with respect to the neutral molecule. Note that H 1 and H 2 represent the alcoholic hydrogen atoms.

		1790 0.0726 0.1258 0.2174 0.0730 0.1457 0.2162 0.0703 0.1433 0.2218 0.0677 0.1448
	H 2	0.1791 0.0726 0.1258 0.2174 0.0740 0.1457 0.2163 0.0702 0.1433 0.2217 0.0677 0.1447
	O	0.1146 0.2682 0.1914 0.1176 0.2670 0.1923 0.1202 0.2684 0.1943 0.1214 0.2642 0.1928
	C	0.1277 0.2126 0.1702 0.1115 0.2050 0.1583 0.1139 0.2086 0.1613 0.1109 0.2199 0.1654
	H 3	0.0428 0.1249 0.0839 0.0471 0.1251 0.0861 0.0475 0.1251 0.0863 0.0489 0.1260 0.0875
	H 4	0.0429 0.1249 0.0839 0.0470 0.1250 0.0860 0.0475 0.1250 0.0862 0.0489 0.1260 0.0875
	H 5	0.0388 0.1951 0.1170 0.0437 0.1907 0.1172 0.0443 0.1911 0.1177 0.0458 0.1887 0.1173

Table 3 . 5 :

 35 Condensed-to-atom Fukui indices for nucleophilic (f + ), electrophilic (f -), and radical attack (f 0 ) for the methoxide evaluated with different XC functionals. In red are highlighted the highest values of the indices for each specific kind of chemical attack: the oxygen atom appears to be the most reactive site for all kind of chemical attacks.

		0186 0.3894 0.2040 0.0399 0.3572 0.2095 0.0151 0.3958 0.2045 0.0210 0.3958 0.2084
	C	0.0004 0.0108 0.0056 0.0014 0.0051 0.0033 0.0001 0.0002 0.0002 -0.0001 0.0008 0.0004
	H 3	0.0046 0.1147 0.0597 0.0115 0.01184 0.0650 0.0040 0.1244 0.0642 0.0058 0.1255 0.0657
	H 4	0.0046 0.1149 0.0598 0.0115 0.1187 0.0651 0.0039 0.1246 0.0643 0.0058 0.1256 0.0971
	H 5	0.0046 0.1151 0.0599 0.0115 0.1188 0.0652 0.0039 0.1248 0.0644 0.0058 0.1259 0.0659

which are a soft acid and an hard base, respectively. This distinction is able to account for the relative production yields. Indeed, since "soft likes soft and hard likes hard"

[START_REF] Pearson | Hard and Soft Acids and Bases[END_REF]

the first reaction path (i.e., that leading to formaldehyde, methane, and water) is more favored than the second one. Of course, as shown in Fig.

3

.6-b, formaldehyde, being an extremely reactive molecular species, can be further employed in the dimethyl ether production. As a corollary, several different calculations within this theorethical context have been performed; all of them carry the same informations of the Fukui analysis (see Appendix B) and more or less obvious quantitative descriptions about the global molecular electrophilic power (e.g., see Table

B

.4 of Appendix B) and the other global and local quantities.

Table B . 3 :

 B3 Table B.2: Chemical potential in eV of the analyzed molecular species. Global chemical softness in eV of the analyzed molecular species.

				Chemical potential µ (eV)
	Molecule		PBE	BLYP	PW91	PZ (LDA)
	CH 3 OH CH 3 OH 2 + CH 3 O -	-5.4883 -9.4597 -2.1888	-5.7502 -10.3616 -2.3385	-5.4943 -10.3862 -2.1901	-5.5388 -10.4289 -2.2416
					Softness S (eV)
	Molecule	PBE	BLYP	PW91	PZ (LDA)
	CH 3 OH	0.1130	0.1217	0.1120	0.1129
	CH 3 OH 2	+	0.0818	0.0773	0.0763	0.0775
	CH 3 O -		0.5389	0.4203	0.5317	0.5141
				Electrophilic Power ω (eV)
	Molecule		PBE	BLYP	PW91	PZ (LDA)
	CH 3 OH		3.4038	4.0232	3.3813	3.4621
	CH 3 OH 2	+		7.3161	8.2952	8.2267	8.4318
	CH 3 O -			2.5818	2.2983	2.5505	2.5834

Table B . 4 :

 B4 Global electrophilic power in eV of the analyzed molecular species. As expected, the electrophilic character of the methanol cation CH 3 OH 2 + is considerably the higher one whereas, on the other hand, the nucleophilic tendency of methoxide is clearly evident by its ω minimum value.Table B.5:Local philicity indices (ω k,α =ωf k,α ) for nucleophilic (ω +,α ), electrophilic (ω -,α ), and radical attack (ω 0,α ) of the atomic sites α of the investigated samples. Note that H 1 and H 2 represent the alcoholic hydrogen atoms of the methanol cation and, at the same time, H 1 is the alcoholic hydrogen of the methanol molecule.

	Local philicity		CH 3 OH			CH 3 OH 2	+	CH 3 O -
	Atom	ω +	ω -	ω 0	ω +	ω -	ω 0	ω +	ω -	ω 0
	H 1	0.0671	0.3118	0.1896	1.3096	0.5311	0.9204
	H 2				1.3096	0.5311	0.9204
	O	0.0333	1.5688	0.8011	0.8384	1.9622	1.4003	0.0480	1.0054	0.5167
	C	0.0238	0.2886	0.1562	0.9343	1.5554	1.2452	0.0010	0.0279	0.0145
	H 3	0.0197	0.5490	0.2846	0.3131	0.9138	0.6138	0.0119	0.2961	0.1541
	H 4	0.0330	0.3002	0.1668	0.3139	0.9138	0.6138	0.0119	0.2966	0.1544
	H 5	0.0197	0.5497	0.2849	0.2839	1.4274	0.8560	0.0119	0.2972	0.0599
	Local hardness	H 1	H 2	O	C	H 3	H 4	H 5
	CH 3 OH	0.5549		2.8506	0.5211	0.9961	0.5404	0.9971
	CH 3 OH 2	+	0.1051	0.1051	2.8393	2.0855	1.3697	1.3697	2.2653
	CH 3 O -				0.5731	0.0160	0.1697	0.1700	0.1703

Table B . 6 :

 B6 Local hardness indices (η α ) for the α atoms of the three investigated molecular species. Note that H 1 and H 2 represent the alcoholic hydrogen atoms of the methanol cation and, at the same time, H 1 is the alcoholic hydrogen of the methanol molecule.

Acknowledgements

Table 3.1: Percent fractions of the molecular species present at the end of each field exploration. The lack of formaldehyde above the field threshold which is able to create it (i.e., 0.55 V/ Å) does not means that it is not present in a remarkable amount during the total simulation. Indeed formaldehyde acts as an intermediate species which contributes to the complexification of the system.

which has been possible to create. However, after 2.5 ps of dynamics performed for each of these two field intensities, only the existence of water, dimethyl ether, and methoxymethanol (in order of decreasing amount) has been possible to observe. Formaldehyde in fact is a sort of intermediate species which mainly and ultimately contributes to the process of polymerization/complexification; it appears the most reactive species among the neutral ones, as it is already known.

Because of the lack of traces of formaldehyde at the end of our numerical experiment (i.e., at 0.75 V/ Å), the mean lifetime of this species has been evaluated at different electric field strengths, as shown in Table 3.2. As a result, it is no possible to find a relationship between these two quantities and what is clear from the analysis of the trajectories is that the local environment (i.e., the solvent) is much more decisive on the reactivity of a given species than the field strength. This latter acts as a catalyst in the sense that it can manipulate the pH of a solution: in all the observed chemical reactions it is evident the key role played not only by the methanol cations/anions but also by the hydronium and hydroxide species obtained as a Field strength (V/ Å) Observed reaction mechanism 0.00

Table 3.6: Observed mechanisms of reaction 3.1 in liquid in absence (metadynamics driven reaction) and in presence of a field strength (0.55 V/ Å) that makes spontaneous the reaction.

Conclusions

In this thesis the molecular dissociation and the proton transfer phenomena have been exstensively investigated in different H-bonded systems.

In particular, analyzing the effect of proton ordering on the transport properties of water ices I h and XI under intense electric fields allowed us to find that the enhanced coupling between the dipole moments and the internal electric field makes it simpler for water molecules to dissociate in the ferroelectric phase than in ice I h . On the other side, restoring the dipole moments distribution after a proton has migrated via the Grotthuss mechanism is less problematic in ice I h than in ice XI, where all molecular dipoles need to be re-tilted in order for another proton to be transferred along a given path. This leads to a higher static conductivity of common ice as compared with that of its proton-ordered counterpart. This "ambivalent" role played by proton (dis)order was an unexpected and so far unreported feature of water dissociation and proton conduction in ices. Concerning this specific aspect, we conclude that orientational order favors, on an energetic basis, the dissociation of water molecules but hinders the proton diffusion which, instead, is promoted and amplified by proton disorder with its associated contribution to the entropy of the system. Moreover, the fact that the dissociation and current thresholds are the same in liquid water and ice I h suggests that the stability of the molecular interactions does not affect the conductivity and that the symmetry of the problem plays a crucial role in the protolysis as well as in the proton transfer phenomena.

123

In order to complete the study on the topological dependencies associated with the proton migration, we have investigated some aspects concerning the dynamics of the proton transfer mechanism in ice I h and ice XI under the action of a static electric field which has been oriented along two different directions in both samples. After obtaining the equation of state, we have calculated the mechanical response through the stiffness tensor. We notice that it does not depend of the field intensity -at least in the investigated regime -when the field is applied along the c-axis of the I h phase. Moreover, also ice XI shows the same trend when the field is oriented along the b-axis. This is no longer true when the field is parallel to the ferroelectric axis. In fact, a monotonic decrease of the elastic constant c 66 has been observed when the field intensity increases, which suggests a softening of the structure.

In addition to the mechanical behavior of the samples, we also have focused on their electrical responses. As far as the ice XI is concerned, as expected from its intrinsic anisotropy, some differences are observed when different directions of the external field are chosen. In particular, the one-stage mechanism, involving at the same time molecular dissociation and proton transfer, observed when the field direction coincides with the ferroelectric axis, disappears when the field is oriented along another axis (e.g., the b-axis), in which case the process develops in two stages.

Moreover, for high field strengths dipole moments tend to align with the field in ice XI more easily than in the I h phase because of the ferroelectric nature exhibited by many of its planes along the b-axis direction. This pseudo-ferroelectricity is enhanced and propagated by fieldinduced defects formation.

A significant finding of this study is the acquisition that the underlying mechanism of molecular dissociation is substantially the same in liquid as well as in solid water, involving a high degree of cooperation of the oxygen sublattice. The dynamics of the molecules is characterized by a compression of the H-bond wires that participate in the dissociation events, independently of the symmetry of the system. Another common feature shared by all the investigated ice systems is the increase of the length of the H-bond chain which immediately follows its compression, a

Appendix A Theoretical background on Fukui

Functions, (Regio)Selectivity, and

Reactivity Descriptors

A.1 Global reactivity descriptors

As it is well-known in DFT the electron density plays the protagonist role since in principle from its knowledge can be possible to determine any physical observable. However, when dealing with the reactivity of an atom or of a molecule, the absolute value of the electron density is not enough by itself. Indeed during a chemical reaction the structural changes or the variations of the external potential are more important in the evaluation of the route that the process will take.

Actually, there exist many global parameters that are somewhat the reflection of the response of the electron density to a given well specified perturbation. The chemical potential µ, electronegativity χ, hardness η, and softness S are defined as follows [START_REF] Parr | Absolute Hardness: Companion Parameter to Absolute Electronegativity[END_REF]:

where of course E, N , and v are the energy, the number of electrons, and the external potential of the molecular system under consideration, respectively. Hence, µ and η represent the slope and curvature of the energy with respect to the electron number variation. Since this latter can only change in a discrete manner, the practical way to face the discontinuity is to apply the finite difference approximation. Moreover, since all the just defined quantities are evaluable at fixed external potential v, a geometry optimization with fixed parameters of the molecule and of its electronic anionic and cationic counterparts is mandatory. This way, once the energies of the neutral molecules and of its "electronically perturbed" analogues are known, it is straightforward to evaluate all the global reactivity descriptors as follows:

where IE and EA are the ionization energy and the electron affinity, respectively. In the frozen orbital method these values can be associated to the highest occupied molecular orbital (HOMO) and to the lowest unoccupied one (LUMO) energies.

Although all these descriptors have been used inter alia to give a theoretical justification to the well-known hard and soft acids and bases (HSAB) principle [START_REF] Pearson | Hard and Soft Acids and Bases[END_REF][START_REF] Pearson | Recent Advances in the Concept of Hard and Soft Acids and Bases[END_REF] and it exists a fundamental link between softness and polarizability [START_REF] Ghanty | Correlation between Hardness, Polarizability, and Size of Atoms, Molecules, and Clusters[END_REF][START_REF] Roy | Correlation of Polarizability, Hardness, and Electronegativity: Polyatomic Molecules[END_REF][START_REF] Ayers | The Physical Basis of the Hard/Soft Acid/Base Principle[END_REF], their use can give only an indication on the reactivity behavior of a global molecular system. demonstrated its usefulness as an intramolecular reactivity descriptor for electrophilic aromatic substitution, addition reactions, and alkylation and acylation of lithium enolates [START_REF] Meneses | A Proposal for a New Local Hardness as Selectivity Index[END_REF].

Finally, we want to stress that still exists an amount of doubtfulness concerning some criteria about the FF. Their importance can be indeed elucidated via several pertubarbation analysis which led to somewhat clouded conclusions about their employment. In particular, it was believed that for hard-hard reactions the most reactive sites were those showing the minima values of the Fukui indices [START_REF] Li | The Fukui Function: a Key Concept Linking Frontier Molecular Orbital Theory and the Hard-Soft-Acid-Base Principle[END_REF], whereas it has been more recently clarified [START_REF] Chattaraj | Chemical Reactivity and Selectivity: Local HSAB Principle versus Frontier Orbital Theory[END_REF] that this minimum condition belongs to a more general situation in which the most preferred site is that which contains the maximum net charge. This sometimes may coincide with the minimum criterion. The soft-soft reactions are orbital controlled and there exists a total consesus about the reliability of a maximum Fukui index criterion which reflects itself in the maximization of the other local reactivity and selectivity descriptor parameters: the most reactive site is that which shows the maxima FF values. Since conceptual (or chemical) DFT has to face with concepts that are quantitative but at same time may be not perfectly rigorous, one must to pay more attention in analyzing each quantity. For this reason a systematic study by exploiting also different levels of DFT approximations should be recommended in some cases.

A.3 Another approach: the Noncovalent Interactions

Though this method, to our best knowledge, has never been used before in conjunction with the FF approach or with the global reactivity descriptors, we propose to use it in order to visualize better the behavior of the electronic density and of its derivatives during the chemical reactions under investigation. After all, the Noncovalent Interactions (NCI) approach shares some similarities with the previously described methods despite the obvious requirement of a single point calculation. Indeed, again the quantum-mechanical electron density n(r) plays the unique and maybe irreplaceable role of protagonist, as almost always occurs in DFT. The

Appendix