Marica Di

Andrea Mele

Frederik Philippi

Ines Adriaan

Layla Luke

Nicolas

Ryan Mirella

Oliver

Laura

Sachini, Tarek Fernando

Stéphane Parola

Carine Michel

Edwige, Damien Marie-Françoise I Tao

J S Freeman

K Goloviznina

H Li

M Saunders

G G Warr

Keywords: polarisation, force field, molecular dynamics, ionic liquids, eutectic solvents. v polarisation, champ de forces, dynamique moléculaire, liquides ioniques, solvants eutectiques. vi

First of all, I would like to express my gratitude to my supervisors, Agilio Padua and Margarida Costa Gomes. Agilio, thank you for your patience and your readiness to discuss my results and to answer my questions any time I needed, for guiding me in the right direction, for motivating me, for respecting my opinion, and for giving me freedom of choice. Margarida, thank you for your valuable suggestions on my research, for introducing me to experimental physical chemistry, for your useful pieces of advice on how to present in front of the public, and for sharing your experience not only in the professional area but in daily life. I appreciate your kindness, your continuous support, and your contribution to my personal growth as a scientist. I was very fortunate to have both of you as supervisors.

My deepest thanks to all the members of the jury for their willingness to evaluate my thesis. To Céline Merlet and to Luka Monticelli, who agreed to come to the ENS Lyon despite a complicated sanitary situation. To Edward Maginn and to Paul Madden, who participated virtually.

This project could not be

Structural formula C 1 C 1 im + 1,3-dimethylimidazolium C 2 C 1 im + 1-ethyl-3-methylimidazolium C 4 C 1 im + 1-butyl-3-methylimidazolium C 6 C 1 im + 1-hexyl-3-methylimidazolium Npmim + 1-methyl-3-neopentyl-imidazolium
Ionic liquids (ILs) are salts formed by large organic ions with delocalized electrostatic charge, conformational flexibility and asymmetric molecular shapes, which hinder crystallisation and lower the melting point to near room temperature [START_REF] Marcus | Room Temperature Ionic Liquids[END_REF][START_REF] Rogers | Ionic Liquids -Solvents of the Future?[END_REF]. They are represented by the imidazolium, pyridinium, pyrrolidinium, quaternary ammonium and phosphonium families with different substituents, which can be combined with a large variety of anions, ranging from small, symmetrical tetrafluoroborate BF - 4 to big, flexible bistriflamide NTf - 2 (Fig. 1.1). Due to their unique characteristics, such as low flammability, negligible vapour pressure, high conductivity, good thermal and electrochemical stability [START_REF] Hallett | Room-Temperature Ionic Liquids: Solvents for Synthesis and Catalysis. 2[END_REF][START_REF] Welton | Room-Temperature Ionic Liquids. Solvents for Synthesis and Catalysis[END_REF], ILs are considered as promising solvation media. Moreover, the properties of ILs can be easily tuned though combining different cations and anions, and functionalizing side-chains, which makes them suitable compounds for effective solvent design. For instance, replacing the chloride anion with triflate (OTf -) in 1-ethyl-3-methylimidazolium ILs leads to a reduction in viscosity from 341 mPa s -1 to 20 mPa s -1 at 323 K [START_REF] Seddon | Viscosity and Density of 1-Alkyl-3-methylimidazolium Ionic Liquids[END_REF], which has a positive impact on conductivity.

Ionic liquids find their application as solvents for reactions, separations, synthesis of materials, as catalysts, as electrolytes in energy-storage devices, as lubricants, optical and magnetic fluids, etc. [START_REF] Ghandi | A Review of Ionic Liquids, Their Limits and Applications[END_REF]. They are known to dissolve a broad range of solutes including carbon dioxide [START_REF] Muldoon | Improving Carbon Dioxide Solubility in Ionic Liquids[END_REF], cellulose [START_REF] Swatloski | Dissolution of Cellose with Ionic Liquids[END_REF] and water insoluble dyes [START_REF] Castner | Intermolecular Dynamics, Interactions, and Solvation in Ionic Liquids[END_REF], making them attractive media for greenhouse gas capturing [START_REF] Lepre | Ionic Liquids Can Enable the Recycling of Fluorinated Greenhouse Gases[END_REF], textile recycling [START_REF] Meksi | A review of progress in the ecological application of ionic liquids in textile processes[END_REF], biomass processing [START_REF] Tan | Ionic Liquids in Biomass Processing[END_REF], and also in pharmaceutical and cosmetic industries [START_REF] Marrucho | Ionic Liquids in Pharmaceutical Applications[END_REF].

ILs show a variety of intra-and intermolecular interactions including electrostatic, van der Waals, π-stacking, H-bonding components, etc. Due to their ionic nature, they are often compared to high-temperature molten salts, and as liquid salts they display charge alternation and intermediate range ordering [START_REF] Sharma | A Brief Guide to the Structure of High-Temperature Molten Salts and Key Aspects Making Them Different from Their Low-Temperature Relatives, the Ionic Liquids[END_REF]. However, the structural patterns may contain different features, for example, microphase separation between polar and non-polar domains in ILs, associated to charged head groups and long alkyl side-chains (Canongia Lopes et al., 2006b). Charge delocalisation and large ion size lead to weaker interactions between ions, attenuating the long range ordering, and possible cation-anion hydrogen bonding results in structural directionality, especially pronounced for ILs in their crystalline phases [START_REF] Dupont | From Molten Salts to Ionic Liquids: A "Nano"Journey[END_REF]. ILs are typically less dense than molten salts because they contain mainly light atoms, with density values in the range 0.9-1.4 g cm -3 for ILs at room temperature [START_REF] Marcus | Room Temperature Ionic Liquids[END_REF], compared to 1.7-1.9 g cm -3 for molten salts at 750-1300 K [START_REF] Redkin | Electrical Conductivity of Molten Electrolytes with a Common Cation[END_REF]. The larger ion size and higher viscosity of ILs, which ranges from tens to thousands of mPa s at 298 K, have an impact on the electrical conductivity, which is orders of magnitude lower than that of hightemperature molten salts. For example, the conductivity of LiCl at 1200 K is 7420 mS cm -1 [START_REF] Redkin | Electrical Conductivity of Molten Electrolytes with a Common Cation[END_REF], while for [C 4 C 1 im][Ntf 2 ] the value is only 4 mS cm -1 [START_REF] Kanakubo | Temperature and Pressure Dependence of the Electrical Conductivity of 1-Butyl-3-methylimidazolium Bis(trifluoromethanesulfonyl)amide[END_REF]. Conductivity of ILs depends on ion paring [START_REF] Kashyap | How Is Charge Transport Different in Ionic Liquids and Electrolyte Solutions?[END_REF], the fraction of ions available for the charge transport (ionicity) has typical values of 0.5-0.8 [START_REF] Macfarlane | On the concept of ionicity in ionic liquids[END_REF][START_REF] Nordness | Ion Dissociation in Ionic Liquids and Ionic Liquid Solutions[END_REF][START_REF] Ueno | Ionicity in ionic liquids: correlation with ionic structure and physicochemical properties[END_REF].

Protic ionic liquids (PILs, Fig. 1.2) are a subclass of ionic liquids formed by a Brønsted acid and a Brønsted base [START_REF] Greaves | Protic Ionic Liquids: Solvents with Tunable Phase Behavior and Physicochemical Properties[END_REF]. Upon mixing, transfer of a proton from acid to base takes place resulting in a purely ionic fluid, whereas the backward transfer is a rare event in the liquid state, provided a sufficient acidity difference. For example, the self-dissociation constant of ethylammonium nitrate (EAN), the most well characterized PIL, equals pK s = 10.0 at 298 K [START_REF] Kanzaki | Acid-Base Property of Ethylammonium Nitrate Ionic Liquid Directly Obtained Using Ion-selective Field Effect Transistor Electrode[END_REF], which indicates a small concentration of nitric acid and ethylamine at equilibrium, of 10 -5 mol l -1 . The distinguishing features of PILs are a non-negligible vapour pressure with the possibility to distil some of them, which becomes feasible due to the reverse proton transfer in the gas phase, upon formation of molecular acid and base species (Greaves et al., 2008b). The presence of extensive hydrogen bonds (Hayes et al., 2013) leads to an increase of cohesive energy, stronger ion packing and enhanced ion correlations between the cations and anions [START_REF] Greaves | Protic Ionic Liquids: Evolving Structure-Property Relationships and Expanding Applications[END_REF]. This explains the higher viscosity and lower ionicity of PILs [START_REF] Yaghini | Effect of water on the transport properties of protic and aprotic imidazolium ionic liquids -an analysis of self-diffusivity, conductivity, and proton exchangemechanism[END_REF] when compared to aprotic ILs, although their transport properties remain interesting for many applications (EAN: η = 33 mPa s, κ = 27 mS cm -1 at 298 K [START_REF] Greaves | Protic Ionic Liquids: Solvents with Tunable Phase Behavior and Physicochemical Properties[END_REF]). Like aprotic ILs (Canongia Lopes et al., 2006b), PILs may have a tendency to nanoscale segregation [START_REF] Atkin | The Smallest Amphiphiles: Nanostructure in Protic Room-Temperature Ionic Liquids with Short Alkyl Groups[END_REF].

PIL are often compared to water because of the similar 3D H-bond network [START_REF] Fumino | Hydrogen Bonding in Protic Ionic Liquids: Reminiscent of Water[END_REF] and of the ability to assist micellisation and to support microemulsions and lyotropic liquid crystal phases (Greaves et al., 2008a[START_REF] Greaves | Protic Ionic Liquids: Evolving Structure-Property Relationships and Expanding Applications[END_REF]. Since these media can conduct protons, PILs are suitable as non-aqueous electrolytes in Li-ion batteries [START_REF] Menne | Protic ionic liquids as electrolytes for lithium-ion batteries[END_REF] and fuel cells [START_REF] Yasuda | Protic ionic liquids: Fuel cell applications[END_REF], and are commonly used as catalysts with the activity strongly correlated with acidic strength (Greaves et al., 2008b). Their broad use also includes liquid and gas chromatography, CO 2 and SO 2 absorption, cellulose dissolution and regeneration, organic and inorganic synthesis (Greaves et al., 2008b[START_REF] Greaves | Protic Ionic Liquids: Evolving Structure-Property Relationships and Expanding Applications[END_REF]. When talking about ionic liquids, one should consider deep eutectic solvents (DES) [START_REF] R E F E R E N C E S Abbott | Novel solvent properties of choline chloride/urea mixtures[END_REF], a related class of compounds. Contrary to ILs, DES are not purely ionic, but mixtures of a quaternary ammonium or phosphonium salt with a molecular compound, a hydrogen bond donor (HBD), the most common examples of which are given in Fig. 1. [START_REF]EtOSO ] and [Aliquat][dca][END_REF]. By definition, eutectic mixtures have lower melting point (T m ) than each of the components. A representative DES is choline chloride plus urea (1 : 2 molar ratio) with a T m = 12 • C, which is notably lower that the melting points of its individual compounds, T m (ChCl) = 302 • C and T m (U) = 133 • C. This lowering of the melting point may arise from a decrease of the cation-anion electrostatic interactions by the HBD, together with steric effects, frustrating crystallisation [START_REF] Hansen | Deep Eutectic Solvents: A Review of Fundamentals and Applications[END_REF]. The majority of DES are known to have relatively high viscosities (> 100 mPa s) at room temperature, which is attributed to an important hydrogen bond network, large ion size and small free volume, being also reflected in relatively low ionic conductivities (< 2 mS cm -1 ) [START_REF] Hansen | Deep Eutectic Solvents: A Review of Fundamentals and Applications[END_REF]Zhang et al., 2012a). The viscosity of the choline chloride urea (1 : 2) mixture reaches 1398 mPa s at 298 K, and conductivity is only 0.4 mS cm -1 [START_REF] Agieienko | Densities, Viscosities, and Electrical Conductivities of Pure Anhydrous Reline and Its Mixtures with Water in the Temperature Range[END_REF].

Ease of preparation from cheap and commercially available chemicals makes DES advantageous systems for industrial usage. Thanks to a good solvation ability, DES are used in electrodeposition of metals, preparation of materials and nanoparticles, gas adsorption, biotransformations and organic synthesis [START_REF] Hansen | Deep Eutectic Solvents: A Review of Fundamentals and Applications[END_REF][START_REF] Liu | Natural Deep Eutectic Solvents: Properties, Applications, and Perspectives[END_REF][START_REF] Paiva | Natural Deep Eutectic Solvents -Solvents for the 21st Century[END_REF]. As multi-component systems, their properties can also be modified by changing the composition. Nevertheless, DES cannot completely replace conventional ILs because of relatively high viscosities caused by hydrogen bonding, high reactivity and instability during electrochemical processes [START_REF] Smith | Deep Eutectic Solvents (DESs) and Their Applications[END_REF]Zhang et al., 2012a), which make them a class of solvents with somewhat different domains of application from aprotic ILs.

rational design of il and des by means of classical md simulations

A huge variety of ILs and DES can be obtained via ion replacement and side-chain modification and cannot be fully characterized by experimental techniques, so theoretical methods are indispensable in their rational design according to the principles of sustainable chemistry [START_REF] Welton | Solvents and sustainable chemistry[END_REF]. Molecular modelling provides access to properties that are difficult or impossible to reach experimentally and contributes to the understanding of structural organisation of the nano-scale level and to the nature of physical and chemical processes [START_REF] Lynden-Bell | Simulations of Ionic Liquids, Solutions, and Surfaces[END_REF]. For example, the domain nanostructure of ILs was initially predicted by all-atom simulations (Canongia Lopes et al., 2006b) and then confirmed by X-ray scattering [START_REF] Triolo | Nanoscale Segregation in Room Temperature Ionic Liquids[END_REF].

Modelling of ionic liquids and related compounds is a challenging task due to the diversity of their interactions, a subtle balance of which determines the behaviour of such complex fluid phases. A choice of a good method capable of representing medium-range ordering (10 1 -10 2 nm) and of capturing structural changes related to ion diffusion (10 1 -10 2 ns) is important for a complete and reliable description of ILs.

Static quantum chemical (QC) calculations and ab initio molecular dynamics (AIMD) are powerful techniques to study local electronic and quasi-chemical effects [START_REF] Kirchner | Multiresolution Calculation of Ionic Liquids[END_REF]. These methods were used to describe hydrogen bonding in ILs and DES, contributed to investigation of the proton transfer [START_REF] Campetella | Unexpected proton mobility in the bulk phase of cholinium-based ionic liquids: new insights from theoretical calculations[END_REF][START_REF] Low | An ab initio Study of the Structure and Energetics of Hydrogen Bonding in Ionic Liquids[END_REF][START_REF] Zhu | Proton Transfer in Phosphoric Acid-Based Protic Ionic Liquids: Effects of the Base[END_REF], and revealed possible carbene formation in [C 2 C 1 im][OAc] [START_REF] Brehm | Proton transfer and polarity changes in ionic liquid-water mixtures: a perspective on hydrogen bonds from ab initio molecular dynamics at the example of 1-ethyl-3-methylimidazolium acetate-water mixtures-Part 1[END_REF][START_REF] Hollóczki | Carbene Formation in Ionic Liquids: Spontaneous, Induced, or Prohibited?[END_REF]. However, such a detailed level of description is reflected in a high computation cost of the methods, and their application is limited to systems composed of only dozens of ion pairs simulated for a few hundreds of picoseconds (Fig. 1.4), this being insufficient to predict of energetic, structural and dynamic properties of the bulk liquids, given the slow dynamics and strong ordering. Much greater time and size scales can be achieved in coarse-grained (CG) simulations, which do not treat all atoms explicitly, but consider a group of atoms in every CG site or bead. For example, the 1-ethyl-3-methylimidazolium cation consists of three beads, a charged head group and two alkyl side-chains, as shown in Fig. 1.5. Despite the simplified representation, the CG force fields can predict nano-scale segregation of bulk ILs (Bhargava et al., 2007b) and describe their behaviour at electrode surface [START_REF] Merlet | New Coarse-Grained Models of Imidazolium Ionic Liquids for Bulk and Interfacial Molecular Simulations[END_REF]. They are indispensable for modelling highly ordered systems, usually with slow dynamics, such as ionic liquid crystals and polymers [START_REF] Wang | Coarse-grained simulations of ionic liquid materials: from monomeric ionic liquids to ionic liquid crystals and polymeric ionic liquids[END_REF], or surface active ILs that can aggregate and form micelles [START_REF] Li | Aggregation Behavior of Trisiloxane-Tailed Surface Active Ionic Liquids in Aqueous Solution: Coarse-Grained Molecular Dynamics Study[END_REF].

Since interaction parameters in CG models lose their atomic character and are obtained though a force-matching procedure, these force fields become system-specific, and so cannot be used together with other models [START_REF] Salanne | Simulations of room temperature ionic liquids: from polarizable to coarsegrained force fields[END_REF]. The majority of CG force fields for ILs cover only the imidazolium family [START_REF] Fajardo | Molecular dynamics simulation of imidazolium CnMIM-BF4 ionic liquids using a coarse grained force-field[END_REF][START_REF] Moradzadeh | Coarse-Grained Force Field for Imidazolium-Based Ionic Liquids[END_REF][START_REF] Wang | Transferable Coarse-Grained Models for Ionic Liquids[END_REF], and their extension to include other ions requires full parametrisation. Transferability can be achieved if a fragment-based approach is adopted as in the MARTINI (Marrink et al., 2007; [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] force field. The MARTINI force field is based on a mapping, according to which three-four heavy atoms plus hydrogens attached are merged into one interaction site with polarity and hydrogen-bonding capability characters. The force field has been recently extended to ionic liquids [START_REF] Vazquez-Salazar | Martini coarse-grained models of imidazolium-based ionic liquids: from nanostructural organization to liquid-liquid extraction[END_REF], but it may not be suitable for DES due to the small size of the interacting species, for which a higher lever of detail may be desirable.

In united-atom force fields, all heavy atoms are treated explicitly, and hydrogen atoms are merged to the atoms they are bounded to. Non-polarisable [START_REF] Liu | A novel united-atom force field for imidazolium-based ionic liquids[END_REF][START_REF] Martínez-Jiménez | United atom model for ionic liquids: UAM-IL[END_REF][START_REF] Urahata | Structure of ionic liquids of 1-alkyl-3-methylimidazolium cations: A systematic computer simulation study[END_REF][START_REF] Zhong | Improved Classical United-Atom Force Field for Imidazolium-Based Ionic Liquids: Tetrafluoroborate, Hexafluorophosphate, Methylsulfate, Trifluoromethylsulfonate, Acetate, Trifluoroacetate, and Bis(trifluoromethylsulfonyl)amide[END_REF] and polarisable models [START_REF] Son | First-Principles United Atom Force Field for the Ionic Liquid BMIM+BF4-: An Alternative to Charge Scaling[END_REF] are available in the literature and provide reasonable description of IL properties, but still lack a good representation of hydrogen bonding, which is particularly important for protic ionic liquids and deep eutectic solvents as this type of interaction determines their unique properties.

All-atom molecular simulation can be considered as a compromise between level of detail and computational cost, attaining size and time scales of 6-10 nm and 10-100 ns, respectively. In the last decades, general non-polarisable force fields, such as the All-Atom Optimized Potentials for Liquid Simulations (OPLS-AA) [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF], the General AM-BER Force Field (AMBER/GAFF) [START_REF] Wang | Development and testing of a general amber force field[END_REF] or the CHARMM General Force Field (CHARMM/CGenFF) [START_REF] Vanommeslaeghe | CHARMM General Force Field: A Force Field for Druglike Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields[END_REF], were used to model ILs and related systems. Individual terms such as partial charges for ions or torsion energy profiles were often reparametrized to fit ab initio calculations, but at that time experimental properties such as densities, viscosities or conductivities were scarce compromising the applicability of these models to ILs.

Design of IL-specific fixed-charge force field for ionic liquids started from the works of Lynden-Bell [START_REF] Hanke | Intermolecular potentials for simulations of liquid imidazolium salts[END_REF] and of Maginn [START_REF] Morrow | Molecular Dynamics Study of the Ionic Liquid 1n-Butyl-3-methylimidazolium Hexafluorophosphate[END_REF], who concentrated their efforts on [C n C 1 im][PF 6 ], and were followed by systematic force fields by Wang [START_REF] Liu | A Refined Force Field for Molecular Simulation of Imidazolium-Based Ionic Liquids[END_REF], Canongia Lopes and Padua (Canongia Lopes et al., 2004a,b, 2006a, 2008, 2012;[START_REF] Shimizu | Molecular Force Field for Ionic Liquids V: Hydroxyethylimidazolium, Dimethoxy-2-Methylimidazolium, and Fluoroalkylimidazolium Cations and Bis(Fluorosulfonyl)Amide, Perfluoroalkanesulfonylamide, and Fluoroalkylfluorophosphate Anions[END_REF], and Acevedo [START_REF] Doherty | Revisiting OPLS Force Field Parameters for Ionic Liquid Simulations[END_REF][START_REF] Sambasivarao | Development of OPLS-AA Force Field Parameters for 68 Unique Ionic Liquids[END_REF]. The CL&P model by Canongia Lopes and Padua, which covers imidazolium, pyridinium, ammonium and phosphonium families with different functional groups and a broad series of anions, became one of the most well-known and widely-used force fields for ILs due to its compatibility with OPLS-AA, robustness and reliability, and numerous efforts to refine it took place during the last decade [START_REF] Chaban | Polarizability versus mobility: atomistic force field for ionic liquids[END_REF][START_REF] Dommert | Refining classical force fields for ionic liquids: theory and application to [MMIM][Cl[END_REF][START_REF] Zhao | A Refined All-Atom Model for the Ionic Liquid 1-n-Butyl 3-Methylimidazolium bis(Trifluoromethylsulfonyl)imide [bmim] [Tf2N][END_REF]. For deep eutectic solvents, a more recent domain than ILs, the representative models are those of Colina [START_REF] Perkins | Molecular Dynamic Simulations and Vibrational Analysis of an Ionic Liquid Analogue[END_REF][START_REF] Perkins | Experimental and Computational Studies of Choline Chloride-Based Deep Eutectic Solvents[END_REF], Aparicio [START_REF] Alcalde | An experimental and theoretical investigation of the physicochemical properties on choline chloride -Lactic acid based natural deep eutectic solvent (NADES)[END_REF][START_REF] García | The impact of charges in force field parameterization for molecular dynamics simulations of deep eutectic solvents[END_REF] and Acevedo [START_REF] Doherty | OPLS Force Field for Choline Chloride-Based Deep Eutectic Solvents[END_REF], and nowadays this field evolves rapidly.

Depending on the assumptions of a model, the net ionic charge in a non-polarisable force field can be set as integer or fractional. The models with unit ionic charges provide a good insight into the structure of ILs in solid and liquid states (Canongia Lopes et al., 2006b), with deviations from experimental densities not exceeding 2 % for the CL&P model [START_REF] Lopes | CL&P: A generic and systematic force field for ionic liquids modeling[END_REF]. Nonetheless, they systematically and dramatically underestimate diffusion coefficients, often by an order of magnitude [START_REF] Tsuzuki | Molecular Dynamics Simulations of Ionic Liquids: Cation and Anion Dependence of Self-Diffusion Coefficients of Ions[END_REF], failing to predict other dynamic proprieties such as viscosity or conductivity [START_REF] Chaban | Polarizability versus mobility: atomistic force field for ionic liquids[END_REF], and also overestimate the heat of vaporisation [START_REF] Santos | Ionic Liquids: First Direct Determination of their Cohesive Energy[END_REF], indicating a misrepresentation of contributions in the interaction potential.

Improvement of the sluggish dynamics can be achieved though scaling the net charge or modification of the van der Waals potential. In the first approach, denoted as scaled-charge models, the ionic charge is reduced to represent charge transfer and/or polarisation effects observed originally for isolated cation-anion pairs in the gas phase [START_REF] Morrow | Molecular Dynamics Study of the Ionic Liquid 1n-Butyl-3-methylimidazolium Hexafluorophosphate[END_REF]. Distinguishing between these two effects for big organic ions is difficult being a subject of open discussion in the literature [START_REF] Salanne | Simulations of room temperature ionic liquids: from polarizable to coarsegrained force fields[END_REF]. Scaling factors of a typical range of 0.5-0.9 are obtained from ab initio calculations charge distributions on a single ion pair or an ionic cluster [START_REF] Schmidt | Ionic Charge Reduction and Atomic Partial Charges from First-Principles Calculations of 1,3-Dimethylimidazolium Chloride[END_REF][START_REF] Zahn | Charge Spreading in Deep Eutectic Solvents[END_REF]Zhang et al., 2012b), from the electronic dielectric constant [START_REF] Dommert | Force Fields for Studying the Structure and Dynamics of Ionic Liquids: A Critical Review of Recent Developments[END_REF][START_REF] Leontyev | Polarizable molecular interactions in condensed phase and their equivalent nonpolarizable models[END_REF][START_REF] Zhao | A Refined All-Atom Model for the Ionic Liquid 1-n-Butyl 3-Methylimidazolium bis(Trifluoromethylsulfonyl)imide [bmim] [Tf2N][END_REF] or chosen empirically [START_REF] Chaban | Polarizability versus mobility: atomistic force field for ionic liquids[END_REF][START_REF] Chaumont | Aqueous Interfaces with Hydrophobic Room-Temperature Ionic Liquids: A Molecular Dynamics Study[END_REF][START_REF] Youngs | Application of Static Charge Transfer within an Ionic-Liquid Force Field and Its Effect on Structure and Dynamics[END_REF]. For example, Balasubramanian (Bhargava et al., 2007a) proposed a universal value of 0.8 for scaling the charges in the CL&P force field, which increased the individual ionic diffusion coefficients of [C 4 C 1 im][PF 6 ] by one order of magnitude, the relative deviations from experiments not exceeding 20 %. Scaledcharge force fields became a popular approach due to their simplicity, good performance, but their justification by charge transfer effects is still superficial. The choice of scaling factors, derived from QC, for a mixture of ILs with a common ion remains ambiguous [START_REF] Bedrov | Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields[END_REF], and a similar problem appears for deep eutectic solvents, where scaling the partial charges of the neutral HBD component [START_REF] Doherty | OPLS Force Field for Choline Chloride-Based Deep Eutectic Solvents[END_REF][START_REF] Sapir | Restructuring a Deep Eutectic Solvent by Water: The Nanostructure of Hydrated Choline Chloride/Urea[END_REF] or not [START_REF] Ferreira | Improved Force Field Model for the Deep Eutectic Solvent Ethaline: Reliable Physicochemical Properties[END_REF][START_REF] Perkins | Molecular Dynamic Simulations and Vibrational Analysis of an Ionic Liquid Analogue[END_REF] depended on author preferences.

The enhanced dynamics in a system with reduced charges can be explained by the weakening of cation-anion electrostatic interactions, associated to the decrease in cohesive energy. This effect can be interpreted as simulating at higher temperature, with consequent fluidification and volume expansion, which explains a 2-7 % density underestimation. In spite of notable changes in cation-anion interactions, the cation-cation and anion-anion spatial correlations are not affected, with the long-range structure being kept (McDaniel et al., 2018b;Yan et al., 2010b). Modifications of short-range electrostatics also introduce large deviations to the mean rotational relaxation times and produce inaccurate dipole moment distributions [START_REF] Schröder | Comparing reduced partial charge models with polarizable simulations of ionic liquids[END_REF].

Furthermore, a disagreement is observed between simulations and experiments when solvation phenomena are studied. Due to the considerable underestimation of cohesive energy, the scaled-charge force fields fail to predict phase separation in [C 4 C 1 im][BF 4 ] -poly(ethylene oxide) mixtures at any temperature [START_REF] Choi | Entropic Mechanism for the Lower Critical Solution Temperature of Poly(ethylene oxide) in a Room Temperature Ionic Liquid[END_REF]. For mixtures of ILs with low dielectric solvents (dichloroethane, acetone, etc.), the model predicts liquid-liquid separation when complete mixing should be observed [START_REF] Mcdaniel | Polarization Effects in Binary [BMIM+][BF4-]/1,2-Dichloroethane, Acetone, Acetonitrile, and Water Electrolytes[END_REF]. In addition, polar gases in ILs undergo aggregation and not dissolution, with calculated free energy and enthalpy of solvation being systematically higher [START_REF] Cui | Influence of Charge Scaling on the Solvation Properties of Ionic Liquid Solutions[END_REF]. Consequently, the scaled-charge force fields can provide reliable results only in a limited number of cases and should be used with care.

The second approach to modify non-polarisable force fields keeps the ionic charges integer, and use empirical adjustment of the van der Waals potential [START_REF] Bedrov | Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields[END_REF]. Ludwig [START_REF] Köddermann | Molecular Dynamic Simulations of Ionic Liquids: A ReliableDescription of Structure, Thermodynamics and Dynamics[END_REF][START_REF] Köddermann | Ionic Liquids: Dissecting the Enthalpies of Vaporization[END_REF] reparametrized the CL&P force field of [C n C 1 im][NTf 2 ], n = 1, 2, 4, 6 and 8 to improve the prediction of heats of vaporisation, self-diffusion coefficients, viscosity and reorientational correlation times resulting in a good agreement with experimental values. Chaumont [START_REF] Chaumont | Is Charge Scaling Really Mandatory when Developing Fixed-Charge Atomistic Force Fields for Deep Eutectic Solvents?[END_REF] proposed a model for ethaline and glyceline DES with the refined Lennard-Jones parameters of a hydroxyl group that was capable to evaluate dynamic properties in a wide temperature range with relative deviations from experiment below 50 %. Even so, this approach lacks transferability and leads to a system-specific force field, which limits its application.

why are polarisable force fields better?

In order to overcome the shortcomings of fixed-charge force fields (both with unit and fractional charge), electronic effects should be incorporated into all-atom models via explicit polarisation terms [START_REF] Bedrov | Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields[END_REF][START_REF] Salanne | Simulations of room temperature ionic liquids: from polarizable to coarsegrained force fields[END_REF]. This is a simplified representation of electronic clouds, as shown at Fig. 1.6, which is significantly less complex when compared to electronic structure methods. It improves the prediction of transport properties while keeping the same order of computational cost.

The first polarisable model for ionic liquids was proposed almost simultaneously with the scaled-charge non-polarisable models. In 2004, [START_REF] Yan | Molecular Dynamics Simulation of Ionic Liquids: The Effect of Electronic Polarizability[END_REF] developed a polarisable force field for [C 2 C 1 im][NO 3 ] based on the point induced dipole approach. In this pioneering study, the author observed a shift of the cation-anion radial distribution function (RDF) maximum to shorter distances, a disordering of anion-anion interaction and a viscosity decrease caused by the explicit induction effects. Later, a detailed analysis of the same system (Yan et al., 2010a,b) confirmed the strengthening of short-range electrostatic interactions and the attenuation of long-range correlations, revealing changes in ion packing and an enhancement of domain aggregation due to the increased screening. Translational motion of individual ions was accel- erated by 2.5-3 times, explained by weaker ion cages. Kato [START_REF] Nakano | A wave-function based approach for polarizable charge model: Systematic comparison of polarization effects on protic, aprotic, and ionic liquids[END_REF] calculated the cage correlation functions for a few imidazolium-based nitrate and chloride ILs, and demonstrated that ions escape from their surrounding cages 10-50 % faster when polarisation effects are present. Rotational motion was also facilitated, but less notably than translational, the corresponding relaxation times decreased only by 20 % [START_REF] Bedrov | Influence of Polarization on Structural, Thermodynamic, and Dynamic Properties of Ionic Liquids Obtained from Molecular Dynamics Simulations[END_REF]. The influence of polarisation on the thermodynamic properties was examined by Bedrov [START_REF] Bedrov | Influence of Polarization on Structural, Thermodynamic, and Dynamic Properties of Ionic Liquids Obtained from Molecular Dynamics Simulations[END_REF], who found a 15-25 % decrease in the heat of vaporisation upon inclusion of the polarisation terms. This effect was related to the strengthening of the cation-anion interactions in the gas phase due to significant induced dipole moments and to the absence of the effective charge screening.

The choice of the force field is always crucial: ILs and DES can be modelled by general polarisable force fields, such as AMBER [START_REF] Wang | Development of Polarizable Models for Molecular Mechanical Calculations II: Induced Dipole Models Significantly Improve Accuracy of Intermolecular Interaction Energies[END_REF] or CHARMM-Drude [START_REF] Lin | Polarizable Force Field for Molecular Ions Based on the Classical Drude Oscillator[END_REF][START_REF] Lin | Further Optimization and Validation of the Classical Drude Polarizable Protein Force Field[END_REF][START_REF] Lopes | Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator[END_REF], but more specific developments, which take into account their ionic nature, are more promising. This involves extension of existing models to the new classes of compounds, on example an adaptation of the AMOEBA [START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF] to ILs in the AMOEBA-IL [START_REF] Starovoytov | Development of AMOEBA Force Field for 1,3-Dimethylimidazolium Based Ionic Liquids[END_REF][START_REF] Vázquez-Montelongo | Current Status of AMOEBA-IL: A Multipolar/Polarizable Force Field for Ionic Liquids[END_REF] force field, or design of new models. In spite of a large number of isolated works for several specific ionic liquids [START_REF] Cavalcante | Polarizability effects on the structure and dynamics of ionic liquids[END_REF][START_REF] Nakano | A wave-function based approach for polarizable charge model: Systematic comparison of polarization effects on protic, aprotic, and ionic liquids[END_REF][START_REF] Schröder | Simulating polarizable molecular ionic liquids with Drude oscillators[END_REF][START_REF] Yan | Molecular Dynamics Simulation of Ionic Liquids: The Effect of Electronic Polarizability[END_REF], there are just a few systematic force fields proposed so far [START_REF] Borodin | Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids[END_REF][START_REF] Mcdaniel | Ab Initio Force Fields for Imidazolium-Based Ionic Liquids[END_REF][START_REF] Vázquez-Montelongo | Current Status of AMOEBA-IL: A Multipolar/Polarizable Force Field for Ionic Liquids[END_REF].

The APPLE&P force field by [START_REF] Borodin | Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids[END_REF], which is based on point induced dipoles, was the first polarisable force field covering all main ionic liquid families. The same methodology had been previously used by the author to derive a force field for molecular solvents and lithium-based electrolytes (Borodin et al., 2006a(Borodin et al., ,b, 2009a)). Parameterized entirely from scratch, it showed a good prediction ability in thermodynamic, structural and transport properties of pure ILs, with deviations not exceeding 3 % in density, 60 % in diffusion coefficients, and 100 % in conductivity and in viscosity. However, an extension of this force field to other classes of compounds might be difficult because of an uncommon exponential-6 potential and the Waldman-Hagler combining rule used to describe the van der Waals interactions. The repulsion-dispersion parameters were fitted to reproduce experimental density and diffusion coefficients, explaining the good performance of this model in evaluating these properties.

Later, Yethiraj and McDaniel [START_REF] Choi | First-Principles, Physically Motivated Force Field for the Ionic Liquid [BMIM][BF4[END_REF][START_REF] Mcdaniel | Ab Initio Force Fields for Imidazolium-Based Ionic Liquids[END_REF]McDaniel et al., , 2018a) ) proposed a first-principle SAPT-FF force field for imidazolium-based ILs, with polarisation represented by Drude induced dipoles. All intermolecular interactions were parametrized from symmetryadapted perturbation theory (SAPT) calculations, and the intramolecular interaction are either taken from the literature for cations and rigid symmetrical anions (BF - 4 ,PF - 6 ), or fitted to ab initio data for flexible anions (NTf - 2 , FSI -). A complex functional form, with multiple interaction terms includes an exponential Born-Mayer repulsive term and a series of C n /r -n , n = 6, 8, 10, 12 powers for the dispersion term [START_REF] Mcdaniel | Physically-Motivated Force Fields from Symmetry-Adapted Perturbation Theory[END_REF], so its usage in common molecular dynamics (MD) codes is not immediate. The prediction ability of this force field is slightly inferior to that of APPLE&P, with larger deviations from experiments being observed in the calculated properties. The SAPT-FF model was extended to polyethylene oxide [START_REF] Son | Conformational and Dynamic Properties of Poly(ethylene oxide) in BMIM+BF4-: A Microsecond Computer Simulation Study Using ab Initio Force Fields[END_REF] to study solvation of polymers in ILs, and recently to a choline chloride plus urea DES [START_REF] Jeong | A Transferable Polarizable Force Field for Urea Crystals and Aqueous Solutions[END_REF][START_REF] Jeong | Deep Eutectic Solvents: Molecular Simulations with a First-Principles Polarizable Force Field[END_REF] being the second polarisable force field for DES in the literature after the model presented in this thesis.

The key disadvantages of these models are the lack of transferability, poor extendability and considerable computational cost of development. We decided to suggest a new polarisable force field for ionic liquids, deep eutectic solvents and electrolytes, which can be derived from existing non-polarisable models without the need of expensive parametrisation and can be thus easily used by the simulation community.

objectives and strategy of the thesis

The aim of this work is to propose a state-of-the-art molecular model for ionic fluids, which accounts explicitly for induction effects and can be used for in silico design of sustainable solvents and reaction media. The model is expected to be general and transferable, covering all major families of ionic liquids and related compounds; to be easily extendable to new systems via a simple and computationally-light procedure; and to demonstrate a good prediction ability.

The organisation of this manuscript is the following:

chapter 2 gives an overview of the theoretical basics of molecular dynamics simulations and explores the ways to incorporate of explicit polarisation effects. The methods of evaluation of transport properties (viscosity and diffusion coefficients) are discussed.

chapter 3 establishes a line of action to develop the CL&Pol polarisable force field for aprotic ionic liquids via the transformation of the existing CL&P fixed-charge model. Explicit polarisation effects are represented using Drude induced dipoles, and a consequent modification of the non-bonded attractive potential in order to avoid double counting of induction effects is justified.

chapter 4 is devoted to the investigation of gas solvation phenomenon in newly synthesized imidazolium-based ionic liquids with Si-functionalized side-chains. We showcase 1.4 objectives and strategy of the thesis the flexibility and robustness of our model when considering novel systems for solvent development.

chapter 5 presents the applicability of the CL&Pol approach to nanomaterials illustrating the ease of extendability of the model to systems described by fixed-charge force fields with functional forms compatible with OPLS-AA. We reveal the role of explicit induction terms in interfacial effects observed for ionic liquids on graphene.

chapter 6 describes the extension of the CL&Pol force field to electrolytes, protic ionic liquids, deep eutectic solvents, and glycols. Modelling these systems is complex due to the presence of small, highly charged ions or strong hydrogen bonds, which cause trajectory instabilities due to the pull exerted on the induced dipoles. We introduce a charge-dipole damping function to smear interactions with these problematic sites and prevent the "polarisation catastrophe".

chapter 7 reports a study of a prototypical system for Li-ion batteries, which consists of the choline chloride -urea DES with lithium chloride salt and small amounts of water. We demonstrate that the CL&Pol is compatible with existing Drude-based polarisable force field and can be straightforwardly mixed without any additional adjustments.

2

P O L A R I S A B L E F O R C E F I E L D S

Context

We will give a short overview of the technical details underlying molecular dynamics simulations with polarisable force fields. Our choice to represent explicit polarisation effects with Drude induced dipoles in the CL&Pol model will be justified. Since the main advantage of polarisable force fields is a correct prediction of transport properties, we will discuss how to compute reliably diffusion coefficients and viscosity, which will be used later for validation of our model.

basics of molecular dynamics simulations

Molecular dynamics (MD) simulations is a computational method based on solving Newton's equations of motion to obtain the phase trajectory of a set of N atoms, from which equilibrium structure and transport properties of the system can be revealed [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Frenkel | Understanding molecular simulation : from algorithms to applications[END_REF]. For an atom i with mass m i , the equations are

v i (t) = d r i (t) dt , (2.1)
a i (t) = d v i (t) dt = d 2 r i (t) dt 2 = F i m i , (2.2)
where r i is the atomic position, v i is the atomic velocity, a i the corresponding acceleration, and F i is the force acting on atom i derived from a potential energy U( r 1 , r 2 , ..., r N ),

F i = - ∂U( r 1 , r 2 , ..., r N ) ∂ r i . (2.
3)

The equations of motion are solved numerically using the finite-difference approach, in which the trajectory is discretized with a time step ∆t. In order to derive an integration algorithm, the atomic coordinates are expressed as a Taylor expansion, 2.4) usually truncated after the ∆t 2 term. Evaluation of higher order terms is not common due to a significant increase in computational cost. The equations of motion are integrated after calculation of the forces, and one of the widely used approaches is the velocity Verlet algorithm,

r i (t + ∆t) = r i (t) + d r i (t) dt ∆t + 1 2! d 2 r i (t) dt 2 ∆t 2 + 1 3! d 3 r i (t) dt 3 ∆t 3 + ..., (
r i (t + ∆t) = r i (t) + v i (t)∆t + F i (t) 2m i ∆t 2 , (2.5) v i (t + ∆t) = v i (t) + F i (t + ∆t) + F i (t) 2m i ∆t. (2.6)
It is a modified version of the classic Verlet algorithm [START_REF] Verlet | Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules[END_REF], with a different way of computing velocities; it can be shown that both methods are equivalent. The method is simple, fast, numerically stable and reversible in time, guaranteeing good energy and momentum conservation.

The equations of motion can be also expressed in the Lagrangian formulation,

- d dt ∂L(q, q) ∂ q + ∂L(q, q) ∂q = 0, (2.7)
where the q(t) are the generalized positions and q(t) the generalized velocities, and L is the Lagrangian defined as a difference between the kinetic and the potential energies,

L(q, q) = U kin ( q) -U pot (q). (2.8)
Contrary to the Newton's equations of motion, this approach operates with energies and not forces (which are vector variables) being mathematically convenient for complex systems. The use of generalized coordinates simplifies the notation and allows to describe a system with numerous constraints in a more systematic way.

A molecular force field is composed of a functional form and of the set of parameters used to evaluate the intra-and intermolecular potential energy of a system of atoms,

U(r 1 , r 2 , ..., r N ) = U bonded (r 1 , r 2 , ..., r N ) + U non-bonded (r 1 , r 2 , ..., r N ).
(2.9)

The bonded energy term comprises bond stretching and angle bending, both commonly described by a harmonic potential, and dihedral, and improper torsions, modelled often with a cosine series,

U bonded (r 1 , r 2 , ..., r N ) = bonds ij k r,ij 2 r ij -r 0,ij 2 + angles ijk k θ,ijk 2 θ ijk -θ 0,ijk 2 + + dihedrals ijkl 4 m=1 V m,ijkl 2 1 + (-1) m cos(mφ ijkl ) , (2.10)
where r 0,ij and θ 0,ijk are the equilibrium bond length and angle, with k r,ij and k θ,ijk corresponding force constants, and V m,ijkl is the dihedral coefficients. The non-bonded contribu-

2.1 basics of molecular dynamics simulations tion includes the electrostatic interaction between the atomic partial charges q, represented by the Coulomb potential, and the van de Waals term set by the 12-6 Lennard-Jones (LJ) potential,

U non-bonded (r 1 , r 2 , ..., r N ) = atoms i atoms j =i U LJ,ij + U Coul,ij , (2.11) U Coul,ij = 1 4π 0 q i q j r ij , (2.12
)

U LJ,ij = 4 ij σ ij r ij 12 - σ ij r ij 6 .
(2.13)

The latter accounts for the London dispersion and Debye induction (through the well-depth ), and Pauli repulsion (through the repulsive atomic diameter σ). The functional form presented above refers to the OPLS-AA [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF] force field, one of the most well-known models for organic molecules.

The force field parameters are derived from high level first principle calculations or through fit to experimental data (x-ray and neutron scattering, infrared, Raman or nuclear magnetic resonance (NMR) spectroscopy). Then, they can be refined to reproduce better macroscopic properties. The OPLS-AA model was validated through the comparison of calculated densities, heats of vaporisation, heat capacities, molar volumes of organic liquids with the experimental results, providing us with a good and reliable tool for modelling liquid state systems.

In our model, geometric-mean mixing rules are used to evaluate the LJ pairwise parameters between atoms of different type,

ij = √ i j , (2.14)
σ ij = √ σ i σ j , (2.15)
Intramolecular interactions of non-bonded type (LJ and Coulomb) between atoms connected by 1 covalent bond (1-2) or 2 (1-3) are excluded, while 1-4 interactions between atoms separated by 3 bonds are reduced by 0.5. Therefore, for the atoms separated by 3 bonds, the Coulomb and the Lennard-Jones potentials account for a part the torsion energy profile, and the dihedral bonded term should be adjusted carefully to represent only the remaining contribution.

In order to assess bulk properties from MD simulations, periodic boundary conditions should be imposed, and a spherical truncation at r cut is applied to reduce the computational cost. Truncation of potentials that do not decay to zero at r cut leads to accumulation of systematic errors in the energy. To overcome this, a long-range correction should be considered, (2.16) where ρ is the average number density. This is valid only for bulk homogeneous fluids, for which one assumes g(r) ≈ 1 at r > r cut , and modelling of heterogeneous systems should be done carefully [START_REF] Jablonka | Applicability of Tail Corrections in the Molecular Simulations of Porous Materials[END_REF]. A typical cut-off radius for the Lennard-Jones potential r cut 2.5σ (for example, σ C = 3.5 Å corresponds to r cut = 8.75 Å), so that the value of the pair potential at this distance represents only 1-2 % of its well depth, leading to a relatively small correction. The upper limit of the cut-off radius should not exceed half of the length of the simulation box according to the minimum image convention. Nevertheless, a long-range correction, similar to Eq. (2.16) is not applicable to the terms that decays slower than r -3 , such as Coulomb potential, for which the long-range part should be properly evaluated. The potential energy of a system of N atoms, in which each charge i interacts with all periodic images of a box with length L, can be written as

U tail = Nρ 2 ∞ r cut 4πr 2 U(r)dr,
U Coul = 1 2 m∈Z 3 N i N i q i q j r ij + mL .
(2.17)

The Ewald method provides an efficient way to compute this sum. Assuming that each point charge is surrounded by a diffuse charge distribution of the opposite sign and with the same magnitude, the long-range electrostatic energy can be computed in reciprocal space using a Fourier series, with a correction for self-interaction applied. The resulting Coulomb potential is expressed as

U Coul = U short-range + U long-range -U self-ewald = = 1 2 N i N j =i q i q j erfc( √ αr ij ) r ij + 1 2V k =0 4π k 2 |ρ( k)| 2 exp(-k 2 /4α) -(α/π) 1/2 N i=1 q 2 i , (2.18)
where erfc(x) is the complementary error function, ρ i (x) the charge distribution represented with a Gaussian of width 2/α, and k = (2π/L) l is the reciprocal vector of a periodic cubic box with length L and volume V, with a triplet of integers l = (l x , l y , l z ). The computational efficiency for large systems can be increased via incorporation of fast Fourier transform (FFT) into the Ewald summation, as done in the particle-particle particlemesh (PPPM) algorithm. The "particle-particle" part corresponds to the short-range interactions evaluated in real space, and the "particle-mesh" part to the long-range interactions in reciprocal space with the charge density interpolated on a grid decreasing the calculation cost of long-range interactions from O(N 2 ) to O(N log N).

adding polarisation effects

Explicit polarisation effects can be incorporated into classical MD simulations through several techniques, namely point induced dipoles, Drude induced dipoles or fluctuating point charges, as schematically represented in Fig. 2.1.

Point induced dipoles

In the point induced dipole model [START_REF] Ren | Consistent treatment of inter-and intramolecular polarization in molecular mechanics calculations[END_REF][START_REF] Ren | Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation[END_REF][START_REF] Rick | Potentials and Algorithms for Incorporating Polarizability in Computer Simulations[END_REF], mathematical dipoles are placed at the centre of atoms. The dipole µ induced on an atomic site i is propor- tional to the electric fields generated by the partial charges E q,i and the induced dipoles E µ,i ,

µ i = αi • ( E q,i + E µ,i ), (2.19)
where αi is atomic polarisability tensor, often approximated by a scalar α i considering isotropic environments of atoms in molecules.

The polarisation energy U pol comprises charge-dipole U qµ , dipole-dipole U µµ and selfpolarisation contributions U self ,

U pol = U qµ + U µµ + U self = - N i µ i • E q,i - N i µ i • E µ,i + 1 2 N i µ i • α-1 i • µ i . (2.20)
Thus, the force field requires only one additional parameter with a clear physical sense, the atomic polarisability. Including high order multipole electrostatics, beyond dipoles, to improve an accuracy of the expansion, is also possible, but long-range corrections become slow affecting the simulation time. Charge transfer effects, on the other hand, cannot be handled with this model [START_REF] Bedrov | Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields[END_REF][START_REF] Rick | Potentials and Algorithms for Incorporating Polarizability in Computer Simulations[END_REF].

Examples of force fields using point induced dipoles include the general AMBER [START_REF] Wang | Development of Polarizable Models for Molecular Mechanical Calculations II: Induced Dipole Models Significantly Improve Accuracy of Intermolecular Interaction Energies[END_REF] and AMOEBA [START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF] models, and more specialized developments for molten salts by Salanne et al. (2012a,b) and, for ionic liquids, APPLE&P by [START_REF] Borodin | Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids[END_REF] and AMOEBA-IL [START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF][START_REF] Vázquez-Montelongo | Current Status of AMOEBA-IL: A Multipolar/Polarizable Force Field for Ionic Liquids[END_REF], which can be used in the AMBER [START_REF] Salomon-Ferrer | An overview of the Amber biomolecular simulation package[END_REF], OpenMM [START_REF] Eastman | OpenMM 7: Rapid development of high performance algorithms for molecular dynamics[END_REF], and Tinker [START_REF] Lagardère | Tinker-HP: a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields[END_REF] MD packages.

Fluctuating point charges

In the fluctuating point charge model [START_REF] Bauer | Recent Applications and Developments of Charge Equilibration Force Fields for Modeling Dynamical Charges in Classical Molecular Dynamics Simulations[END_REF][START_REF] Rappe | Charge Equilibration for Molecular Dynamics Simulations[END_REF][START_REF] Rick | Dynamical fluctuating charge force fields: Application to liquid water[END_REF], the values of the atomic charges are allowed to change during the simulation mimicking the polarisation effects. The energy required to create a charge q on an atom i can be written as a Taylor series expansion,

U(q i ) = U 0 i + χ 0 i q i + 1 2 J ii q 2 i , (2.21)
For an atom with a charge up to ±1e, the Taylor coefficients χ i and J ii /2 are the electronegativity and the hardness, which can be evaluated from the ionisation potential, IP, and the electron affinity, EA,

χ i = (IP i + EA i )/2, (2.22) J ii = IP i -EA i . (2.23)
The total potential energy of a system of atoms is described by

U(q 1 , q 2 , ..., q N ) = N i U 0 i + χ 0 i q i + 1 2 J ii q 2 i + N i N j>i J ij (r ij )q i q j . (2.24)
Atomic partial charges are computed through energy minimisation, which amounts to an electronegativity equalisation,

χ i = ∂U ∂q i = χ 0 i + J ii q i + N j =i J ij (r ij )q j , (2.25) ∂U ∂q 1 = ∂U ∂q 2 = ... = ∂U ∂q i = ... = ∂U ∂q N . (2.26)
The total charge of the system is constrained as

Q tot = N i=1 q i . (2.27)
Charge transfer can be controlled using specific techniques [START_REF] Nistor | A Generalization of the Charge Equilibration Method for Nonmetallic Materials[END_REF][START_REF] Verstraelen | The Electronegativity Equalization Method and the Split Charge Equilibration Applied to Organic Systems: Parametrization, Validation, and Comparison[END_REF], but, on the other hand, it is difficult to exclude it completely. In addition, out-ofplane polarisation is absent, and the force field cannot be easily mixed with the non-polarisable models.

Because of these disadvantages, the fluctuating charges technique did not gain success in the modelling community being mainly represented by the CHARMM-FQ [START_REF] Patel | CHARMM fluctuating charge force field for proteins: I parameterization and application to bulk organic liquid simulations[END_REF] force field, and only a few unrelated studies are found in the literature for ionic liquids [START_REF] Cavalcante | Polarizability effects on the structure and dynamics of ionic liquids[END_REF][START_REF] Wu | The electronegativity equalization method fused with molecular mechanics: a fluctuating charge and flexible body potential function for [Emim][Gly] ionic liquids[END_REF].

Drude induced dipoles

The Drude induced-dipole method [START_REF] Drude | The theory of optics[END_REF]Lamoureux et al., 2003a,b;[START_REF] Maaren | Molecular Dynamics Simulations of Water with Novel Shell-Model Potentials[END_REF] is nearly equivalent to point induced dipoles [START_REF] Schmollngruber | Comparing induced point-dipoles and Drude oscillators[END_REF], but is based on a different way to materialize the induced dipoles. Similarly to point dipoles, Drude dipoles model requires only atomic polarisability to set up the force field, but an additional particle per polarisable atom is introduced. Both methods predict similar trends in structure and dynamics, some minor deviations being observed for highly polarisable atoms [START_REF] Schmollngruber | Comparing induced point-dipoles and Drude oscillators[END_REF].

Among the available Drude models, SEM-Drude [START_REF] Sharma | SEM-Drude Model for the Accurate and Efficient Simulation of MgCl2-KCl Mixtures in the Condensed Phase[END_REF][START_REF] Sharma | A Brief Guide to the Structure of High-Temperature Molten Salts and Key Aspects Making Them Different from Their Low-Temperature Relatives, the Ionic Liquids[END_REF] for molten salts, SAPT-FF [START_REF] Jeong | A Transferable Polarizable Force Field for Urea Crystals and Aqueous Solutions[END_REF][START_REF] Mcdaniel | Ab Initio Force Fields for Imidazolium-Based Ionic Liquids[END_REF] for ionic liquids and the general CHARMM-Drude [START_REF] Lin | Polarizable Force Field for Molecular Ions Based on the Classical Drude Oscillator[END_REF][START_REF] Lin | Further Optimization and Validation of the Classical Drude Polarizable Protein Force Field[END_REF][START_REF] Lopes | Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator[END_REF] force fields should be mentioned. The Drude method is implemented in the widely used codes, including CHARMM (Lamoureux et al., 2003b), NAMD [START_REF] Jiang | Highperformance scalable molecular dynamics simulations of a polarizable force field based on classical Drude oscillators in NAMD[END_REF], GROMACS [START_REF] Lemkul | Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model[END_REF], LAMMPS [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF] and OpenMM [START_REF] Huang | Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation, validation, and benchmarks[END_REF]. Its suitability to MD and versability make this technique the most obvious choice for our CL&Pol force field. [START_REF] Drude | The theory of optics[END_REF] is formed by two point charges of opposite sign, a positively-charged Drude core (DC), at the centre the atomic site r, and a negatively-charged Drude particle (DP), connected to its core by a harmonic spring with an equilibrium distance d of zero (Fig. 2.2). The pair of charges gives rise to an induced dipole under an external electrostatic field E,

µ i = q D,i d i = q 2 D,i k D E, (2.28)
with magnitude determined by the Drude charges, q D , and by the force constant of the spring, k D , which are related through the polarisability

α i = q 2 D,i k D . ( 2

.29)

A displacement of a DP from its DC should not exceed a typical value of | d| = 0.10-0.15 Å being sufficiently small compared to any interatomic distance, allowing the µ to be considered as a point dipole. The Coulomb interaction of the DC-DP pairs with atomic partial charges and induced dipoles then can be evaluated via point charges only,

U Coul = N i N j>i (q i -q D,i )(q j -q D,j ) | r i -r j | + (q i -q D,i )q D,j | r i -r j -d j | + q D,i (q j -q D,j ) | r i + d i -r j | + + q D,i q D,j | r i + d i -r j -d j | , (2.30)
where q i is total atomic charge, and (q i -q D,i ) is the DC charge. The self-polarisation energy corresponds to the harmonic potential of the Drude springs,

U self = N i 1 2 k D d i • d i . (2.31)
We opt, as in the CHARMM polarisable force field by Lamoureux and Roux (2003b), to set the force constant of the DC-DP harmonic bond to k D = 4184 kJ mol -1 Å -2 and to obtain the Drude charge from the atomic polarisability. Alternatively, the Drude charge can be fixed (for example, q D = -1.0e), and the force constant can be evaluated from Eq. (2.29), as proposed by [START_REF] Schröder | Simulating polarizable molecular ionic liquids with Drude oscillators[END_REF], but this may lead to a significant elongation of a DC-DP bond, invalidating the point-dipole approximation. Atomic polarisabilities derived specifically for ionic liquids, from first principle calculations, are available from the literature (Heid et al., 2018a,c) and are listed in Tab 2.1. We opted to use those values instead of general polarisabiltites for organic compounds [START_REF] Kang | Additivity of atomic static polarizabilities and dispersion coefficients[END_REF][START_REF] Miller | Additivity methods in molecular polarizability[END_REF] because the former are more accurate and take into consideration the ionic nature of our systems. Alternatively, molecular (ion pair) polarisability can be obtained from the experimental refractive index via the Lorenz-Lorenz equation, and then decomposed into atomic contributions using regression methods [START_REF] Bica | Polarisabilities of alkylimidazolium ionic liquids[END_REF]. This requires an additional correction of ±1.87 Å 3 spread over the atoms of each ion, to account for the negative or positive unit charge [START_REF] Bernardes | Additive polarizabilities in ionic liquids[END_REF]. Hydrogen atoms are treated as non-polarisable in our model because of low polarizability (low electron density) and small atomic mass, compared to the mass of the DP. Their polar-izability is merged onto the heavy atoms to which they are bound (Lamoureux et al., 2003b;[START_REF] Schmollngruber | Comparing induced point-dipoles and Drude oscillators[END_REF].

Dynamics of Drude particles

In a MD trajectory, the degrees of freedom (DOF) associated with Drude dipoles can be relaxed iteratively to their equilibrium positions at each time step, following the self-consistent field (SCF) condition,

∂U ∂ d i = ∂U self ∂ d i + ∂U elec ∂ d i = 0. (2.32)
So, the equilibrium force condition is expressed via .33) But this self-consistent relaxed-Drude method is rather slow, with the number of required interactions depending on a convergence criterion δ,

k D d i = q D,i E. ( 2
∂U ∂ d i < δ, (2.34)
the weakening of which introduces uncontrolled noise propagation over the simulation. When the minimisation is not completely converged, a systematic drag force appears, having a strong impact on stability of trajectory and calculated properties [START_REF] Heid | Polarizable molecular dynamics simulations of ionic liquids: Influence of temperature control[END_REF]Lamoureux et al., 2003b;[START_REF] Rupakheti | Statistical mechanics of polarizable force fields based on classical Drude oscillators with dynamical propagation by the dual-thermostat extended Lagrangian[END_REF].

A faster alternative is to integrate the equations of motion of the DP by assigning a mass to them, which is subtracted from the respective DC. The aim of this extended Lagrangian method is to generate a trajectory as close as possible to that of the relaxed Drudes. This implies the use of a dual thermostat (Lamoureux et al., 2003b), in which the relative motion of the DP with respect to their DC is kept at a very low temperature [START_REF] Sprik | Computer simulation of the dynamics of induced polarization fluctuations in water[END_REF] (typically of 1 K) to prevent kinetic energy draining into the DC-DP DOF, while the remaining DOF are thermostated at the desired temperature.

In the dual Nosé-Hoover thermostat (NH), implemented in the main MD packages [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF][START_REF] Huang | Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation, validation, and benchmarks[END_REF][START_REF] Jiang | Highperformance scalable molecular dynamics simulations of a polarizable force field based on classical Drude oscillators in NAMD[END_REF][START_REF] Lemkul | Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model[END_REF], the first thermostat keeps the DCs and non-polarisable atoms at the temperature T ,

m i ∂ 2 R i ∂t 2 = F i -m i ∂ R i ∂t ∂η ∂t , (2.35)
where m i is the total mass of a DC-DP pair, R i = (m i r i + m D d i )/m i is the vector defining the centre of mass (COM) of the DC-DP pairs and F i = -∂U/∂ R i is the force acting on the COM.

The second thermostat couples the DP at the temperature T D in reduced coordinates,

m D,i ∂ 2 d i ∂t 2 = F D,i -m D,i ∂ d i ∂t ∂η D ∂t , (2.36)
where m D,i = m D (1 -m D /m i ) denotes the reduced mass, F D,i = -∂U/∂ d i is the force acting on the reduced mass. The η and η D are specific energies of the heat bath, (2.38) where N f = 3N -N cons is the number of DOF associated with the motion of the non-polarisable atoms or the COM of the DC-DP pairs, and with the number of constrains imposed by SHAKE for X-H bonds taken into account. N D,f = 3N D is the number of DOF associated with the DP motion,

Q ∂ 2 η ∂t 2 = N i m i ∂ R i ∂t 2 -N f k B T, (2.37) Q D ∂ 2 η D ∂t 2 = N i m D,i ∂ d i ∂t 2 -N D,f k B T D ,
Q = N f k B T τ 2 and Q D = N D,f k B T D τ 2 D
are the inertia factors of the NH thermostat related to the characteristic times τ and τ D . The forces can be expressed through the absolute coordinates of the DCs ( r i ) and the DP ( r D,i = r i + d i ),

F i = - ∂U ∂ r i - ∂U ∂ r D,i , (2.39)
F D,i = -1 - m D m i ∂U ∂ r D,i + m D m i ∂U ∂ r i . (2.40)
The integration of the equations of motion is split into a one-step velocity-Verlet propagation and a multi-step Nosé-Hoover propagation, which helps to avoid systematic errors in η and η D .

The choice of Drude mass and force constant is crucial, because it affects the integration time step. The values of m D = 0.4 u and k D = 4184 kJ mol -1 Å -2 are optimal, allowing for a time step of 1 fs ensuring good energy conservation and simulation stability (Lamoureux et al., 2003b). We set the thermostat coupling constants to τ = 0.2 ps and τ D = 0.05 ps following recent recommendations [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF].

The same dual thermostatting scheme can be applied to the Langevin thermostat, which was implemented in LAMMPS [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF], OpenMM [START_REF] Huang | Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation, validation, and benchmarks[END_REF], and NAMD [START_REF] Jiang | Highperformance scalable molecular dynamics simulations of a polarizable force field based on classical Drude oscillators in NAMD[END_REF],

m i ∂ 2 R i ∂t 2 = F i -m i γ ∂ R i ∂t + f i , (2.41) m D,i ∂ 2 d i ∂t 2 = F D,i -m D,i γ D ∂ d i ∂t + f D,i , (2.42)
where γ and γ D are the friction coefficients (inverse to the coupling times τ and τ D ). The random forces, 2.44) follow the fluctuation-dissipation theorem, with Gaussian random veriable by R(t) and R D (t).

f i = 2γk B T m i 1/2 R(t), (2.43) f D,i = 2γ D k B T D m D,i 1/2 R D (t), (
Contrary to the Nosé-Hoover thermostat, the dual-Langevin thermostat is not recommended for temperature control in polarisable simulations. It can describe correctly the structural properties, but fails when predicting dynamics (diffusion coefficients, conductivities) revealing a strong dependence on the thermostat coupling constants. Moreover, too weak coupling leads to a poor temperature control of the DP motion, and T D can reach 20-70 K [START_REF] Heid | Polarizable molecular dynamics simulations of ionic liquids: Influence of temperature control[END_REF][START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF].

Recently, Yethiraj [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF] has shown that the dual Nosé-Hoover thermostat suffers from significant heat flow from the system to the cold Drude DOF due to overlap of the power spectrum of the Drude dipoles with high-frequency intramolecular modes. This leads to incorrect partitioning of the kinetic energy, the COM translational temperature increases (up to 100 K!) to compensate for the too low temperature of the high-frequency intramolecular modes. This artefact brings systematic deviations in the properties of highly flexible molecules or ions (namely, in ILs and DES) overestimating their diffusion coefficients and underestimating their density. To correct this bias, a separate thermostat should be imposed on the translational motion of the molecular centres of mass. In this triple temperature-grouped Nosé-Hoover thermostat (tgNH) [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF], Eq. (2.37) is split into COM and internal DOF, .47) where M j is the molecular mass, 3N mol are the total translational DOF, and N f -3N mol is the remaining DOFs, the inertia factors equal

Q COM ∂ 2 η COM ∂t 2 = N mol j M j ∂ R COM,j ∂t 2 -3N mol k B T, (2.45)
Q int ∂ 2 η int ∂t 2 = N i m i ∂ R int,i ∂t 2 -(N f -3N mol )k B T, (2.46) ∂ R int,i ∂t = ∂ R i ∂t - ∂ R COM,j ∂t , ( 2
Q COM = 3N mol k B T τ 2 and Q int = (N f -3N mol )k B T τ 2 .
Eq. (2.35), which represents the motion of the ith atom in a jth molecule, is transformed into .48) Thus, the total real particle temperature (T real ), the COM translation temperature (T COM ), 2.50) and the temperature of intermolecular DOF [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF] should be equal to the desired temperature T . The Drude temperature (T Drude ) Thus, we insist that only the triple temperature-grouped Nosé-Hoover thermostat, which was implemented in the OpenMM [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF] and LAMMPS (Goloviznina et al., 2021a) codes, can provide a correct and accurate description of dynamic and structural quantities. A comparative study of the dual NH and the triple tgNH thermostats will be presented in the next chapters.

m i ∂ 2 R i ∂t 2 = F i -m i ∂ R COM,j ∂t • ∂η COM ∂t + ∂ R int,i ∂t • ∂η int ∂t . ( 2
T real = N i m i (∂ R i /∂t) 2 N f k B , (2.49)
T COM = N mol j M(∂ R COM,j /∂t) 2 3N mol k B , (
T Drude = N i m D,i (∂ d i /∂t) 2 3N D k B , ( 2

Damping short-range interactions of Drude dipoles

When two induced dipoles located in the same or different molecules are found too close, the point charge approximation is no longer valid to describe the electrostatics correctly [START_REF] Birge | Calculation of molecular polarizabilities using an anisotropic atom point dipole interaction model which includes the effect of electron repulsion[END_REF][START_REF] Thole | Molecular polarizabilities calculated with a modified dipole interaction[END_REF] leading to instabilities in the simulation, such as collapsing of a DP into a neighbouring DC. This is especially pronounced for the 1-2, 1-3 and 1-4 intramolecular dipoledipole interactions, since they are not excluded (or scaled), unlike charge-charge or van der Waals terms [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF]. In order to represent smearing of electron clouds and prevent divergence of electrostatic potential, damping of interactions between induced dipoles should be introduced at short range.

The damping function should be finite or approach zero at r → 0, and decay fast enough at large distance to recover the unscreened Coulomb law. Several linear and exponential damping functions were examined in the literature [START_REF] Duijnen | Molecular and Atomic Polarizabilities: Thole's Model Revisited[END_REF][START_REF] Taylor | The effect of Thole functions on the simulation of ionic liquids with point induced dipoles at various densities[END_REF][START_REF] Thole | Molecular polarizabilities calculated with a modified dipole interaction[END_REF]Wang et al., 2012a), among which the functional form

S ij (r ij ) = 1 -1 + s ij r ij 2 exp(-s ij r ij ), (2.52
) (2.53) shown in Fig. 2.4, became the most commonly used due to its good flexibility. Its damping range is controlled by the atomic polarizabilities α, and a parameter a ij . The latter parameter is often considered as universal, namely in APPLE&P [START_REF] Borodin | Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids[END_REF], AMOEBA-IL [START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF][START_REF] Vázquez-Montelongo | Current Status of AMOEBA-IL: A Multipolar/Polarizable Force Field for Ionic Liquids[END_REF], and SAPT-FF [START_REF] Mcdaniel | Ab Initio Force Fields for Imidazolium-Based Ionic Liquids[END_REF] etc., but sometimes is fitted during the force field development procedure, like in a recent CHARMM model [START_REF] Lopes | Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator[END_REF]. We adopted a universal value a ij = 2.6 following [START_REF] Noskov | Molecular dynamics study of hydration in ethanol-water mixtures using a polarizable force field[END_REF].

s ij = a ij (α ij ) 1/3 = (a i + a j )/2 [(α i α j ) 1/2 ] 1/3 ,

Performance and system setup

MD simulations of ILs are usually performed on systems consisting of 100-1000 ionic pairs rarely exceeding 30 000 atoms. We set a cut-off radius to 12 Å, which forces the minimal box length of 24 Å according to the minimum image convention. The upper system size is limited by the memory of the machine and the speed of the execution of the MD code, the most notable contribution to which rises from force calculation [START_REF] Allen | Computer simulation of liquids[END_REF], in particular long-range electrostatics. For simulation of bulk liquids, we considered 300 ion pairs to be a good choice for most systems studied herein with reasonable amount of computational resources per time unit simulated.

A system of 300 ion pairs of [C 4 C 1 im][NTf 2 ] consists of 12000 atoms forming a box with equilibrium an length of 53 Å. In LAMMPS [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF], this can be modelled with a performance of 17.8 ns/day on a single Intel 6142 Sky Lake node with 32 central processing unit (CPU) cores, using a fixed-charge force field. Introducing explicit polarisation adds 7500 DP, and the productivity drops to 3.2 ns/day. The average computational cost of 1 ns of polarisable simulation is about 240 CPU hours. In OpenMM [START_REF] Eastman | OpenMM 7: Rapid development of high performance algorithms for molecular dynamics[END_REF], which runs on graphics processing units (GPU), the productivity can reach 35 ns/day with the dual Nosé-Hoover thermostat, and almost 70 ns/day with the dual Langevin thermostat on a single Nvidia Tesla V100 GPU. The differences between the thermostats in LAMMPS are not significant, and, unfortunately, GPU acceleration is not available for polarisable simulations at this time.

We considered simulation times of 2 ns for equilibration and 10-20 ns for production sufficient for good sampling of the phase space of bulk liquids with viscosities lower than 50 mPa s -1 . The equilibration step was controlled through analysis of energetic (total, potential, kinetic energy, etc.) and macroscopic (temperature, density, etc.) properties, their averages becoming independent in time once equilibration is achieved. An example is provided in Fig. 2.5.

dynamic properties as a criterion of reliability of the force field

As it was already mentioned, one of the main advantages of polarisable force fields is their ability to describe better the dynamics of an ionic fluid [START_REF] Bedrov | Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields[END_REF][START_REF] Salanne | Simulations of room temperature ionic liquids: from polarizable to coarsegrained force fields[END_REF], without need of specific adjustments as in scaled-charge models (Bhargava et al., 2007a). Therefore, it is very important for us to establish a clear and well-defined protocol of calculation of transport properties. For the validation of the force field, we have chosen diffusion coefficients and viscosity, experimental data sets of which are present in the literature for ILs and DES, and can be measured in our group for newly synthetized systems.

Diffusion coefficients from Einstein's relation

We calculated diffusion coefficients from mean-squared displacements using Einstein's relation, (2.54) according to the recommended procedure [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF]. A linear fit is performed in the diffusive regime, typically in the range of 2-8 ns for a 10 ns trajectory (Fig. 2.6) with short times less than 2 ns excluded because of possible cage regime. The Yeh-Hummer correction [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF][START_REF] Yeh | System-Size Dependence of Diffusion Coefficients and Viscosities from Molecular Dynamics Simulations with Periodic Boundary Conditions[END_REF]) is used to account for finite-size effects,

D = lim t→∞ 1 6 d dt (r(t) -r(0)) 2 ,
D 0 = D + 2.8373k B T 6πηL , (2.55)
where k B is the Boltzmann constant, L is the length of a cubic box, and η is the viscosity of the system at temperature T .

Viscosity from equilibrium and non-equilibrium MD

Viscosity is an easily accessible physical property for fluids via several experimental techniques, [START_REF] Diogo | Viscosity Measurements on Ionic Liquids: A Cautionary Tale[END_REF] but its evaluation from simulations is not straightforward [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF]. Several approaches for viscosity calculation based on equilibrium and nonequilibrium MD methods are discussed below.

Green-Kubo relation

For polarisable MD simulations, which are about 5 times slower than fixed-charge ones, the equilibrium methods are appealing since the same MD trajectory can be used for calculation of viscosity, density, diffusion coefficients and structure analysis. Nevertheless, equilibrium methods are limited to low-viscous systems of 50 mPa s and less, these values being typical of a few ionic liquids (containing DCA -, TCB -and NTf - 2 anions) at room temperature. In order to validate our force field for ILs with higher viscosities, we increased the simulation temperature should by up to 50-80 • C to fluidify the liquids.

The shear viscosity is evaluated from the pressure tensor using a Green-Kubo relation [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF] ,

η = V kT ∞ 0 p xy (t)p xy (0) dt, (2.56)
which can be averaged over the off-diagonal components p xy , p yz and p zx to improve precision. Theoretically, the autocorrelation function (ACF) should approach zero in the long time limit, with the integral converged to a constant value but, in practice, the integrand becomes noisy at long times resulting in inconsistent viscosity values [START_REF] Zhang | Reliable Viscosity Calculation from Equilibrium Molecular Dynamics Simulations: A Time Decomposition Method[END_REF], as illustrated by Fig. The simplest solution is to limit the integration time by a point when the ACF crosses zero value (a so-called zero time) [START_REF] Zhang | A molecular dynamics study of bulk and shear viscosity of liquid iron using embedded-atom potential[END_REF] at about 100-200 ps, but this cut-off is a result of statistical fluctuations, and ignoring the long-time behaviour of the ACF causes underestimation of viscosity [START_REF] Mouas | Diffusion and viscosity of liquid tin: Green-Kubo relationship-based calculations from molecular dynamics simulations[END_REF]. For example, the predicted viscosity of [C 2 C 1 im][DCA] at 303 K is (8.7 ± 0.6) mPa s, which is below the experimental value of 13.9 mPa s [START_REF] Freire | Thermophysical Characterization of Ionic Liquids Able To Dissolve Biomass[END_REF]. model fit to autocorrelation function.

In order to overcome large fluctuations and eliminate the noise at long times, a model with two exponential functions can be fitted to the ACF [START_REF] Fanourgakis | Determining the Bulk Viscosity of Rigid Water Models[END_REF][START_REF] Guo | Viscosity and stress autocorrelation function in supercooled water: a molecular dynamics study[END_REF],

S f ACF (t) S f ACF (0) = (1 -C) cos(ωt) exp((-t/τ f ) β f ) + C exp((-t/τ s ) β s ), (2.57)
where C, ω, τ f , τ s , β f , β s are the fitted parameters, and C determines the weight of slow and fast relaxation processes. This has been tested on the different water models [START_REF] Fanourgakis | Determining the Bulk Viscosity of Rigid Water Models[END_REF], and provided a reasonable fit, but some deviations were observed at the extremum of the ACF. Contrary to water, the fast oscillation pattern of an IL is more complex (Fig. 2.7), and Eq. ( 2.57) is unable to produce an adequate fit.

An alternative is to use a hybrid approach, where the raw simulation ACF is integrated up to a certain time (a so-called "switch time") and beyond it is replaced with a fitted ACF [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF]. Several exponential function are proposed in the literature to fit the slow relaxation [START_REF] Basconi | Effects of Temperature Control Algorithms on Transport Properties and Kinetics in Molecular Dynamics Simulations[END_REF][START_REF] Fernández | Shear viscosity and thermal conductivity of quadrupolar real fluids from molecular simulation[END_REF][START_REF] Spoel | A systematic study of water models for molecular simulation: Derivation of water models optimized for use with a reaction field[END_REF],

f 0 (t) = exp(-t/τ), (2.58) f 1 (t) = a • exp(-t/τ), (2.59) f 2 (t) = a • exp -t β , (2.60) f 3 (t) = exp (-t/τ) β = exp -t β /α , (2.61) f 4 (t) = a • exp (-t/τ) β = a • exp -t β /α , (2.62)
where τ β is replaced by α to simplify the notation. The functions f 0 (t) and f 1 (t), in which the power depends linearly on time, are not flexible enough to fit the tail of the ACF of an IL, as shown in the Fig. 2.9), where f 3 (t) predicts a viscosity of 52.0 mPa s, significantly overestimating the experimental value of 17.4 mPa s. The corresponding viscosity given by the f 2 (t) is 16.0 mPa s justifying the choice of this fitting function for our work.

f 3 (t) is more marked in [C 2 C 1 im][BF 4 ] (Fig.
We based the choice of the switch time on maximizing the R 2 ,

R 2 = 1 - N i=0 y i -f i (t i ) 2 N i=0 y i -y 2 , (2.63)
as a function of the starting time of a fit (Fig. 2.10), and the determination of a well-defined maximum is accelerated by the usage of the Golden-section search. After having found the switch time, the fit is performed until the end of the simulation, which forces the fitting

S ACF t / ps p xy f 2 (t) f 3 (t) Figure 2.9 -Fit of the ACF tail of p xy pressure component of [C 2 C 1 im][BF 4 ] at 323 K in the range 1.2-
10 000 ps with the f 2 (t) and f 3 (t) exponential functions. Both graphs represent the same data sets with the logarithmic (left) and linear (right) scales of the x-axis.

function to approach zero at longer times, and sufficiently long trajectories of 10 ns or more are required. The ACF points are not equidistant in time naturally lowering the contribution of noisy long-time data to the fit. Integration is performed numerically using the trapezoidal rule. This procedure yielded satisfactory results, and the viscosities of a series of ILs with the associated standard deviation will be reported in the next chapters. model fit to running integral.

Instead of fitting the ACF tail, Maginn [START_REF] Zhang | Reliable Viscosity Calculation from Equilibrium Molecular Dynamics Simulations: A Time Decomposition Method[END_REF], recommends fitting the running integral

η(t) = V kT τ 0 p αβ (τ)p αβ (0) dτ, τ ∈ [0 : t] (2.64)
with a double-exponential function [START_REF] Rey-Castro | Transport Properties of the Ionic Liquid 1-Ethyl-3-Methylimidazolium Chloride from Equilibrium Molecular Dynamics Simulation. The Effect of Temperature[END_REF])

η(t) = Aατ 1 (1 -exp(-t/τ 1 )) + A(1 -α)τ 2 (1 -exp(-t/τ 2 )), (2.65)
where A, α, τ 1 , τ 2 are fitted parameters. Then, the viscosity value is calculated at t → ∞ according to the relation

η ∞ = Aατ 1 + A(1 -α)τ 2 .
(2.66)

An alternative method is proposed by Borodin et al. (2009b),

η(t) = η ∞ (1 -exp(-t/τ s ) β s ), (2.67)
the functional formula of which we found not flexible enough to fit our data. The running integrals of three pressure components of [C 2 C 1 im][DCA] taken from a single trajectory are shown in Fig 2.11. Due to the strong noise effects at long times, none of them reached a plateau. The fits of the data up to a cut-off time (t cut ) of 1300 ps, which is close to the value used in the literature for an IL with a similar viscosity [START_REF] Zhang | Reliable Viscosity Calculation from Equilibrium Molecular Dynamics Simulations: A Time Decomposition Method[END_REF], resulted in an average value (9 ± 3) mPa s. [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF] recommend to assign the t cut to the value when the uncertainty equals 40 % of viscosity, which corresponds to 1700 ps and η = (10 ± 4) mPa s. Herein, we considered this standard deviation to be too large. Nevertheless, the fitting function cannot always provide us with a reliable prediction, and the long-time noise in the data can result in unrealistic behaviour of the fitting. In the example, shown in Fig 2 .12, the p xz component is constantly increasing with time, and the corresponding fit function does not reach the plateau either. The fitting coefficients are listed in Appendix (Tab. A.2). One could propose to constrain the fitting parameter τ 2 (and τ 1 ), which actually improves the prediction giving η = (17.5 ± 1.7) mPa s in average. But this solution is not robust, the choice of the upper τ limit is not straightforward and the τ variables are artificially enforced to take the constrained values.

Maginn recommends [START_REF] Zhang | Reliable Viscosity Calculation from Equilibrium Molecular Dynamics Simulations: A Time Decomposition Method[END_REF] averaging over at least 30-40 independent simulations to suppress the noise effect. But this appears as a huge task in computational cost for polarisable simulations. Here we generated five independent trajectories, of 10 ns each, shown in Fig 2.12, where two of the components reached a plateau after 500 ps whereas the third was still affected by long-time noise. The viscosity is identical to the value evaluated from a single trajectory, η = (17.5 ± 1.1) mPa s, with a slightly improved standard deviation as presented in Tab. A.2. The use of a larger number of trajectories does not seems justified in our context of resource use with respect to the quality of the results. Thus, we do not recommend use of fit of the running integral for viscosity calculations from polarisable simulations, due to the significant computational resources required to eliminate the noise.

Periodic perturbation method

Viscosity can be calculated from non-equilibrium MD simulations using the periodic perturbation method [START_REF] Hess | Determining the shear viscosity of model liquids from molecular dynamics simulations[END_REF], by applying a cosine-shaped acceleration with amplitude A along the x-axis, with periodicity along the z-axis, a x (z) = A cos(kz), (2.68)

k = 2π/l z , (2.69)
which generates a velocity profile in the x direction

v x (z) = V(1 -exp(-t/τ r ) cos(kz)), (2.70) V = A ρ ηk 2 , (2.71) τ r = ρ/(ηk 2 ), (2.72)
as illustrated in Fig. 2.13. The amplitude of the velocity profile, V, can be evaluated at each time step, and the inverse viscosity averaged over the production run,

1 η = V A 4π 2 l 2 z ρ , (2.73)
where ρ is the density and l z is the height of the box. The velocity profile generated is excluded from the kinetic energy and is not coupled with a thermostat to obtain a correct viscosity value. The shear rate

s max = max z ∂v x (z) ∂z = A ρ ηk (2.74)
should not be too high, in order to keep the simulation close to equilibrium. The standard deviation in viscosity can be computed via

σ η = 2 s max k B T η t a V , (2.75)
where t a is the total simulation time and V the system volume.

We evaluated the viscosity with acceleration amplitudes in the range 0.005-0.1 × 10 -5 Å fs -2 . As discussed previously by [START_REF] Hess | Determining the shear viscosity of model liquids from molecular dynamics simulations[END_REF], too large A values lead to high shear rates, resulting in viscosity underestimation. Too small A values, on the other hand, affect the amplitude of the velocity profile and introduce noise in the calculated data, as reported in Appendix (Tab. A.3) and shown in Fig 2 .14. One could suggest to extrapolate the viscosity to zero acceleration amplitude [START_REF] Zhao | On the accuracy of predicting shear viscosity of molecular liquids using the periodic perturbation method[END_REF], but the calculation at a series of A values has a high computational cost when a polarisable force field is used. Thus, we propose to evaluate viscosity at one value of acceleration amplitude, 0.02 × 10 -5 Å fs -2 , found to be optimal for the systems studied.

The periodic perturbation method yields a better prediction with a smaller standard deviation than the fit of the ACF, and thus we considered to be the best option for viscosity calculation. However, in LAMMPS its usage is limited to flexible molecules (Goloviznina et al., 2021a), since rigid bodies are handled by a specific thermostat. Consequently, the viscosity of systems containing, for example, rigid SWM4-NDP water [START_REF] Lamoureux | A polarizable model of water for molecular dynamics simulations of biomolecules[END_REF]) cannot be computed with this approach.

To sum up

The CL&Pol polarisable force field will be based on Drude induced dipoles, with shortrange dipole-dipole interactions smeared using the Thole damping function. The performance of the dual and temperature-grouped Nosé-Hoover thermostats, used for temperature control in polarisable simulations, will be analysed in the next chapters. We propose to compute diffusion coefficients from mean-squared displacements using Einstein's relation. We recommend that viscosity be evaluated using the periodic perturbation method (from non-equilibrium molecular dynamics) or from the pressure tensor using a Green-Kubo relation with the smoothed tail of the ACF (from equilibrium molecular dynamics), if the viscosity is sufficiently low or the former method is not accessible.

3 T R A N S F E R A B L E , P O L A R I S A B L E F O R C E F I E L D F O R A P R O T I C I O N I C L I Q U I D S

Context

Here we will focus on the development of the CL&Pol polarisable force field for aprotic ionic liquids, which will be based on the existing CL&P fixed-charge model. The transformation will require modification of the Lennard-Jones parameters to compensate for new explicit polarisation terms, introduced via Drude induced dipoles.

from a fixed-charge to a polarisable force field

The development of a polarisable force field from scratch is a huge task, undertaken by some groups for ionic liquids [START_REF] Borodin | Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids[END_REF][START_REF] Mcdaniel | Ab Initio Force Fields for Imidazolium-Based Ionic Liquids[END_REF][START_REF] Starovoytov | Development of AMOEBA Force Field for 1,3-Dimethylimidazolium Based Ionic Liquids[END_REF][START_REF] Vázquez-Montelongo | Current Status of AMOEBA-IL: A Multipolar/Polarizable Force Field for Ionic Liquids[END_REF], but transferability remains an issue. Alternatively, one can attempt to add polarisation terms to a fixed-charge force field, which implies the modification of the non-bonded attractive energies (represented by Lennard-Jones or equivalent potentials describing van der Waals interactions) in order to remove a double counting of the induction effects [START_REF] Padua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF], which are already included implicitly in the fixed-charge model. The non-bonded potential, then, should account only for London dispersion, the polarisation being represented explicitly by the induced dipoles.

Several strategies of force field adjustment are proposed in the literature, based either on the full reparametrisation of non-bonded attractive terms [START_REF] Dočkal | Molecular Force Field Development for Aqueous Electrolytes: 2. Polarizable Models Incorporating Crystalline Chemical Potential and Their Accurate Simulations of Halite, Hydrohalite, Aqueous Solutions of NaCl, and Solubility[END_REF][START_REF] Lemkul | An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications[END_REF][START_REF] Mondal | Quantitative Prediction of Physical Properties of Imidazolium Based Room Temperature Ionic Liquids through Determination of Condensed Phase Site Charges: A Refined Force Field[END_REF][START_REF] Uhlig | First-Principles Parametrization of Polarizable Coarse-Grained Force Fields for Ionic Liquids[END_REF]Wang et al., 2012b) or on the fine tuning of existing models [START_REF] Becker | Investigating polarization effects of CO2 adsorption in MgMOF-74[END_REF][START_REF] Bedrov | Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields[END_REF][START_REF] Padua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF].

Holm [START_REF] Uhlig | First-Principles Parametrization of Polarizable Coarse-Grained Force Fields for Ionic Liquids[END_REF] re-evaluated the C 6 molecular attractive term though the relation

C 6 = 3 π ∞ 0 α(iν)α(iν) dν, (3.1)
where α(iν) is the isotropic dynamic polarisability at the frequency ν. Being proportional to the squared fraction of the atomic volume V i to the total molecular volume V tot [START_REF] Gould | How polarizabilities and C6 coefficients actually vary with atomic volume[END_REF], the individual atomic contributions (i) were isolated,

C i 6 = V i V tot 2 C 6 , (3.2)
and the Lennard-Jones parameters were consequently extracted, C i 6 = 4 i σ 6 i .

Vlugt [START_REF] Becker | Investigating polarization effects of CO2 adsorption in MgMOF-74[END_REF] scaled the well-depth, i , of an original non-polarisable Lennard-Jones potential by an empirical factor evaluated from the atomic polarisability α i ,

k i = (1 + λ) -(α i /α max ) (1 + λ) -λ(α i /α max ) , (3.3)
where α max is the largest atomic polarisability in the system and λ is an empirical factor between 0 and 1 that describes the relative importance of dispersive interactions in the potential. The λ parameter is fine-tuned thought comparison of simulation results to experimental data, leading to 0.3 [START_REF] Zeman | A coarse-grained polarizable force field for the ionic liquid 1-butyl-3-methylimidazolium hexafluorophosphate[END_REF], 0.4 (Heid et al., 2018b), and 0.7 [START_REF] Becker | Investigating polarization effects of CO2 adsorption in MgMOF-74[END_REF] values in the literature. The well-depth i of the atom with the highest polarisability is scaled by k i = λ, and all other atoms by k i varying from λ to 1. Eq. ( 3.3) assures that the i of a non-polarisable atom is not modified (α i = 0 Å 3 ). In addition to scaling the well-depth, Holm [START_REF] Zeman | A coarse-grained polarizable force field for the ionic liquid 1-butyl-3-methylimidazolium hexafluorophosphate[END_REF] suggested to apply a uniform correction of 0.9963 to atomic diameters (σ i ) to match the experimental density in his coarse-grained FF for

[C 4 C 1 im][PF 6 ].
Aiming to generate a transferable force field, [START_REF] Padua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF] proposed to scale the LJ terms (well-depth ij ) in a consistent and systematic manner though evaluation of the individual dispersion and induction components using symmetry-adapted perturbation theory (SAPT) [START_REF] Hohenstein | Wavefunction methods for noncovalent interactions[END_REF][START_REF] Jeziorski | Perturbation theory approach to intermolecular potential energy surfaces of van der Waals complexes[END_REF][START_REF] Parker | Levels of symmetry adapted perturbation theory (SAPT). I. Efficiency and performance for interaction energies[END_REF]. In this QM approach, the energy of a dimer is described as a perturbation to a "non-interacting" Hamiltonian, represented by a sum of the Hamiltonians of two monomers. Contrary to supermolecular methods, SAPT computes the interaction energy directly, without the need to add or subtract quantaties, being free of the basis set superposition error (BSSE). Moreover, meaningful components can be isolated from the interaction energy, i.e. electrostatics, exchange-repulsion, induction and dispersion, as is required for our force field development procedure.

Depending on the level of approximation, several classes of SAPT are defined. In the simpliest case, SAPT0, the monomers are treated at Hartree-Fock (HF) level, and second order perturbation is applied only to the dispersion term (equivalent to MP2),

E SAPT0 = E HF + E (20) disp + E (20) exch-disp disp = = E (10) elst elst + E (10) exch exch + E (20) ind,r + E (20) exch-ind,r + δE (2) HF ind + E (20) disp + E (20) exch-disp disp , (3.4)
where in E (AB) X,(r) , X denotes the contribution to the interactions at A intermolecular and B intramolecular orders of perturbation, and the subscript r denotes the contribution due to the orbital relaxation (response) effects; E HF is the Hartree-Fock energy and δE (2) HF is the difference between HF and SAPT interaction energies, which mainly represents polarisation beyond the second order E exch-disp terms, to improve the representation of hydrogen bonded complexes. This method is relatively inexpensive computationally and provides a quick qualitative analysis. We use it for potential energy scans. The accuracy can be increased by a factor of two (and consequently the computational cost by a factor of 50) by using higher level SAPT2+, which considers second order terms for all intra-and intermonomer components,

E SAPT2+ = E SAPT0 + E (12) elst,r elst + E (11) exch + E (12) exch exch + E (22) ind + E (22) exch-ind ind + + E (21) disp + E (22) disp disp , (3.5)
Here the intramonomer corrections to the dispersion energy are equivalent to the MP4 method. Due to its performance and relatively low mean absolute errors, we use this method to evaluate dispersion and induction components at the potential energy minimum, from which scaling factors are derived. Therefore, SAPT allows access to physically motivated scaling factors, but its high computational cost can be a main drawback, making addition of new compounds to the force field laborious.

The aims of this chapter are to extend this strategy and to devise a scheme of general applicability to transform existing fixed-charge force fields into polarisable versions without the need for expensive first-principles calculations, thus enabling the simulation of ILs with a much better prediction of both equilibrium and transport properties.

fragment approach

We report here a set of SAPT calculations to resolve the different contributions to interaction energies on 46 dimers composed of ions, representative of different families of ionic liquids, and also of neutral molecules. Our aim is to generate a sufficiently large data set to develop a predictive scheme. By parameterising separate fragments and not entire ions [START_REF] Lopes | CL&P: A generic and systematic force field for ionic liquids modeling[END_REF], we can achieve the necessary transferability to describe broad families of ions without the need of a specific parameterisation effort for each individual compound. The structural formulae of the ion families considered here are depicted in Fig. 3.1. The fragment approach followed [START_REF] Padua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF] allows to treat a smaller set of building blocks. For example, most cations in ionic liquids are formed by a charged head-group and one or several side-chains with various possible lengths. So, 1-ethyl-3-methylimidazolium

(C 2 C 1 im + ), 1-ethylpyridinium (C 2 Py + ), 1,1-dimethylpirrolidinium (C 1 C 1 pyrr + )
, and tetrame-thylammonium (N + 1111 ) cations were considered as representative head-groups for imidazolium, pyridinium, pyrrolidinium and quaternary ammonium ionic liquids, respectively. Butane (C 4 H 10 ) was used as the representative fragment for alkyl side-chains, and from its atom types, chains of different lengths can be created. From previous work in force field development for ionic liquids [START_REF] Lopes | CL&P: A generic and systematic force field for ionic liquids modeling[END_REF] we learnt that the influence of a charged group along an alkyl side-chain extends up to two bonds along the chain, so the C 2 and H 1 atoms connected to the first carbon C 1 require specific parameters, especially partial charges. Further along an alkyl side-chain, parameters can be considered identical to those of n-alkanes (Fig. 3.2).

Concerning anions, tetrafluoroborate, dicyanamide (DCA -) and bis(trifluoromethanesulfonyl)imide (NTf -

2 ), bis(fluorosulfonyl)imide (FSI -), trifluoromethanesulfonate or triflate (OTf -), and acetate (OAc -) were treated as entire ions. General alkylsulfonates can in turn be built by combining the parameters of the sulfonate head-group with those of an alkyl side-chain (Canongia [START_REF] Lopes | Molecular Force Field for Ionic Liquids IV: Trialkylimidazolium and Alkoxycarbonyl-Imidazolium Cations; Alkylsulfonate and Alkylsulfate Anions[END_REF]. Here, p-toluenesulfonate or tosylate (TsO -), was split into toluene (Tol) and methanesulfonate or mesylate (MsO -), fragments (Fig. 3.

2).

The fragment set considered also includes neutral polar molecules, such as acetonitrile (AN) and dimethyl sulfoxide (DMSO), the slightly polar dimethoxyethane (DME), which is a monomer unit of polyethylene oxide, non-polar hexane (C 6 H 14 ) and aromatic benzene (Bz). We performed SAPT calculations of the potential energy of interaction, for dimers composed of selected fragments of ionic or molecular species, from which we obtained a sufficiently large set of values for the dispersion and induction energies. Then we proceeded to generate MD trajectories of ionic liquids and their mixtures with molecular compounds to calculate condensed-phase equilibrium and transport properties. These are compared to experimental values, in order to assess the soundness of our strategy through the predictive ability of the polarisable force field.

computational methods and force field

Geometries of isolated molecules, ions and subsequently dimers were optimised using dispersion-corrected density functional theory [START_REF] Grimme | A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[END_REF] at the B97-D3/cc-pVDZ level. The potential energy curves for each dimer were calculated using sSAPT0/jaDZ at a series of the distances between the monomers, while keeping their geometries fixed. At the distance corresponding to the potential energy minimum, energies were calculated at the SAPT2+/aDZ level to obtain more accurate values for the dispersion and induction energies. The jaDZ and aDZ correspond to jun-cc-pVDZ and aug-cc-pVDZ basis sets, augmented versions of the correlation-consistent polarised valence double-zeta (cc-pVDZ) basis set. The jun-cc-pVDZ (or aug-cc-pVDZ') is a truncated aug-cc-pVDZ basis set with all diffuse functions removed from hydrogen atoms and the diffuse d-functions removed from non-hydrogen atoms. These basis sets were chosen for different SAPT levels following the recomendation of Sherrill [START_REF] Hohenstein | Density fitting and Cholesky decomposition approximations in symmetry-adapted perturbation theory: Implementation and application to probe the nature of π-πinteractions in linear acenes[END_REF][START_REF] Hohenstein | Wavefunction methods for noncovalent interactions[END_REF][START_REF] Parker | Levels of symmetry adapted perturbation theory (SAPT). I. Efficiency and performance for interaction energies[END_REF], who considered sSAPT0/jaDZ and SAPT2+/aDZ as "bronze" and "silver" standards. SAPT calculations were performed using the Psi4 code [START_REF] Turney | Psi4: an open-source ab initio electronic structure program[END_REF]. Atomic partial charges and dipole moments of molecules and fragments were obtained, on optimised geometries with Gaussian [START_REF] Frisch | Gaussian~16 Revision B.01. Gaussian Inc[END_REF], using the CHelpG [START_REF] Breneman | Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis[END_REF]) method based on MP2/cc-pVTZ(-f) densities. The dipole moment of charged fragments was calculated using coordinates with origin on the centre of mass (the "standard orientation" of Gaussian).

Molecular dynamics simulations of cubic boxes containing 300 ion pairs for pure ionic liquids, or 300(1 -x) ion pairs and 300x molecules for mixtures with x solvent fraction of molecular compound, were performed with LAMMPS [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF]. Initial configurations were generated using the fftool (Padua, 2021b) and packmol [START_REF] Martínez | PACKMOL: A package for building initial configurations for molecular dynamics simulations[END_REF] utilities. A cutoff of 12 Å was considered for the Lennard-Jones potential, with tail corrections for energy and pressure. The particle-particle particle-mesh method (PPPM) was used to evaluate electrostatic energies with an accuracy of 1 × 10 -5 . Bonds terminating in hydrogen atoms were constrained using the SHAKE algorithm. The timestep was of 1 fs. The systems were equilibrated for 2 ns in the NpT ensemble, following which 10 ns production runs were performed in the NpT and NVT ensembles using dual Nosé-Hoover thermostat and barostat. Different temperatures, in the range 298-353 K, were chosen for different systems according to the availability of experimental data for comparison and to the viscosity of the liquids. Pressure was kept at 1 bar in all runs.

The fixed-charge force field that serves as basis for the present development is the CL&P force field (Canongia Lopes et al., 2004a[START_REF] Lopes | CL&P: A generic and systematic force field for ionic liquids modeling[END_REF]Padua, 2021a), with revised Lennard-Jones parameters for fluorinated sulfonylimide anions [START_REF] Gouveia | Ionic liquids with anions based on fluorosulfonyl derivatives: from asymmetrical substitutions to a consistent force field model[END_REF] that yield more accurate liquid densities with respect to previous work [START_REF] Padua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF], in particular for NTf - 2 -based ILs. Molecular compounds were represented by the OPLS-AA force field [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF].

Adding explicit polarisation requires knowledge of atomic polarisabilities, which we took from the recent work of Schröder (Heid et al., 2018c), as explained in Ch. 2 (section 2.3). All heavy atoms were considered as polarisable, while the polarisabilities of hydrogen atoms were merged onto the polarisability of the heavy atoms to which they are bonded. The mass of Drude particles was set at m D = 0.4 a u and the force constant of the harmonic spring between Drude cores and Drude particles set at k D = 4184 kJ mol -1 Å -2 . The partial charges of the DPs were calculated from the polarisabilities according to Eq. (2.29) (Lamoureux et al., 2003b). The total charge of the DC plus the DP is equal to the initial charge of the (nonpolarisable) atom. Thole damping functions, Eq. ( 2.52), [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF][START_REF] Thole | Molecular polarizabilities calculated with a modified dipole interaction[END_REF] were employed to reduce at short range the electrostatic interactions between induced dipoles, in order to avoid excessive correlation between neighbouring Drude particles. Universal value of the parameter a = 2.6 was chosen according to the functional form given in the literature [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF][START_REF] Noskov | Molecular dynamics study of hydration in ethanol-water mixtures using a polarizable force field[END_REF]. The relative motion of DP with respect to their DC was regulated at an equivalent temperature of 1 K using a specific thermostat [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF]Lamoureux et al., 2003b). Input files for LAMMPS with the polarisable force field were prepared using the polarizer and scaleLJ tools [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF][START_REF] Di Pietro | Perturbing structure and dynamics in hydrated reline with lithium chloride[END_REF], which convert non-polarisable input files to polarisable versions. Simulations with Drude induced dipoles can be enable in LAMMPS by activating the DRUDE package [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF]. Dynamic properties, namely diffusion coefficients and viscosities, were evaluated from equilibrium trajectories using the Einstein, Eq. ( 2.54), and Green-Kubo, Eq. ( 2.56), relations, following the recommendations of a recent review [START_REF] Maginn | Best Practices for Computing Transport Properties 1. Self-Diffusivity and Viscosity from Equilibrium Molecular Dynamics [Article v1.0[END_REF] as explained in Ch. 2 above (section 2.4). Diffusion coefficients were calculated from mean-squared displacements using Einstein's relation yielding converged values within the duration of the trajectories. However, the evaluation of shear viscosity using the Green-Kubo relation was more difficult to converge, because the ACF in the integrand becomes noisy at long times and the viscosity values from integrating the three off-diagonal components, p xy , p yz and p zx become inconsistent. We smoothed the tail of the ACF with an exponential decay function. This procedure yielded satisfactory results, and the viscosities we report are the average of the three off-diagonal components, with the associated standard deviation.

evaluation of induction and dispersion

SAPT calculations of the repulsive, electrostatic, induction and dispersion terms of the potential energy of interaction were performed for 46 dimers composed of charged and neutral fragments: 13 cation-anion, 9 cation-neutral, 14 anion-neutral and 10 neutral-neutral fragment dimers. Two examples, a cation-anion pair We evaluated a k ij factor, which is the ratio between the dispersion contribution and the sum of dispersion and induction,

C 2 C 1 im + • • • BF - 4 and a neutral molecule-ion C 4 H 10 • • • BF - 4 ,
k ij = E disp E disp + E ind . (3.6)
This is the factor by which the LJ attractive energy in the non-polarisable force field should be scaled, in order to retain only dispersive terms and exclude implicit polarisation contributions from the LJ potential, since induction will be explicitly represented by the Drude induced dipoles. Therefore, in the polarisable force field, the LJ parameters will be obtained by scaling by a factor k ij those of the CL&P or OPLS-AA models. The scaling factors are evaluated per fragment, so for the ions and molecules that will be considered in applications, different scaling factors may be applied between atom groups within the ions or molecules (Fig. 3.4).

We considered only dimers that can be found at distances characteristic of first solvation or coordination shells, i.e., cation-anion, ion-neutral and neutral-neutral, so we did not perform calculations between same-charge fragments whose interactions are dominated by the repulsive electrostatic terms and are usually found at large separations. The scaling factors from SAPT are derived from gas-phase dimers, so we sought to evaluate the relative dispersion and induction contributions in condensed-phase simulations of ionic liquids and their mixtures with molecular compounds. We performed MD simulations of 11 ionic liquids based on C 2 C 1 im + , C 4 C 1 im + , C 6 C 1 im + and C 4 C 1 pyrr + cations, and BF - 4 , DCA -and NTf - 2 anions. The present data sets on SAPT calculations on dimers and MD simulations on liquid systems considerably expand the previous report [START_REF] Padua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF]) on 8 dimers and 3 ionic liquids.

E c = E IL tot /N -E + tot -E - tot , E LJ = E IL LJ /N -E + LJ -E - LJ , E elst = E IL elst /N -E + elst -E - elst , (3.7)
where the E IL denotes the system total, van der Waals or electrostatic energy averaged over the MD trajectory; and the E + and E -correspond to the total, van der Waals or electrostatic energy of a single cation or anion obtained from additional 10 ns NVT runs at the same temperature in a simulation box with 30 Å side. The induction contribution in the liquid state can be evaluated as the change in electrostatic energy due to introducing Drude particles, taking into account the self-energy of the Drude oscillators (which amounts to the potential energy stored in the DC-DP harmonic bonds),

E ind ≈ E elst -E elst (FixQ) -E self . (3.8)
The fraction of the dispersion energy in the non-bonded interactions in condensed phase can be evaluated through K = E LJ /( E LJ + E ind ), which can be compared with the effective scaling factors from SAPT, averaged per atom according to the fragments involved, k ij , and presented in Tab. 3.2. In general, the K values from the liquid-state simulations were found to be close to the k ij , with differences below 0.1, which means that the SAPT calculations on isolated dimers are good predictors of the energy decomposition in condensed phase. An exception is found for [RC 1 im][DCA] ionic liquids, for which the induction contribution in the liquid phase appears significantly larger than in the gas-phase dimers. We investigated the source of this discrepancy and found that in the optimised dimer geometry, the dicyanamide anion is placed above the imidazolium ring [START_REF] Kempter | The role of hydrogen atoms in interactions involving imidazolium-based ionic liquids[END_REF], in order to minimize the electrostatic energy of the pair, as shown in Fig. 3.5. In the liquid phase, several configurations are possible [START_REF] Weber | Complex Structural and Dynamical Interplay of Cyano-Based Ionic Liquids[END_REF], including anions sitting close to the plane of the imidazolium ring (Fig. 3.6) with hydrogen bonds being formed between the H atoms of the ring and the terminal N atoms of dicyanamide (Fig. 3.7). Thus, we interpret the discrepancy in the values of the k ij scaling factors as a result of the geometry of the isolated dimer, which does not represent the diversity of liquid phase cation-anion structures. Because our aim is to develop a general methodology, we opted to not fine-tune this step of geometry optimisation of the dimers and instead assess the validity of this choice by the ability of the new model to predict equilibrium and transport properties of ionic liquid systems. Using the k ij scaling factors from SAPT calculations on fragment dimers provides thus a means to adapt an existing force field into a polarisable one, as was shown in a previous feasibility report [START_REF] Padua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF]. Predictions of equilibrium and transport properties of ionic liquids and mixtures with molecular compounds will be presented below. The evaluation of scaling factors using SAPT2+ is nevertheless computationally demanding. For example, a single-point calculation of C 2 Py + • • • NTf 2 dimer at SAPT2+/aDZ level requires 75 GB of memory and takes 4.5 days on 16 processors.

E c E LJ E elst E self E ind k ij K [C 4 C 1 im][BF 4 ] 343 K FixQ -

predicting dispersion scaling factors

The computational cost of SAPT calculations led us to propose a general predictive scheme to obtain the k ij scaling factors from simple, readily calculated atomic or fragment proper- ties. The scaling factor, Eq. (3.6), can be expressed through the ratio between induction and dispersion energies,

k ij = 1 + E ind E disp -1 , (3.9)
and these two terms can be evaluated from induction (Debye) and dispersion (London) forces.

The induction energy between two molecular or ionic fragments i and j, averaged over orientations, is given by [START_REF] Hirschfelder | Molecular Theory of Gases and Liquids[END_REF]) 3.10) where the α are the polarisabilities of the fragments, Q their net charges and µ the dipole moments (we calculated the dipole moment of charged fragments using coordinates with origin on the centre of mass). The distance r ij is measured between the fragments. The first term in Eq. (3.10) is expected to dominate for charged fragments and the second for neutral but polar fragments. The dominant term in the dispersion energy between two monomers is 3.11) and it depends on polarisabilities α and ionisation energies I. The ratio between the induction and the dispersion energies takes thus the form

E ind = - Q 2 i α j + Q 2 j α i 2r 4 ij - µ 2 i α j + µ 2 j α i r 6 ij -. . . , (
E disp = - 3 2
I i I j I i + I j α i α j r 6 ij + . . . , (
E ind E disp = c 0 r 2 ij Q 2 i α j + Q 2 j α i α i α j + c 1 µ 2 i α j + µ 2 j α i α i α j , (3.12)
where we consider c 0 and c 1 as coefficients to fit against the SAPT calculations. The distance r ij is considered here to be the equilibrium distance of each dimer. Thus, the k ij factor is evaluated through the relation

k ij = 1 + c 0 r 2 ij Q 2 i α j + Q 2 j α i α i α j + c 1 µ 2 i α j + µ 2 j α i α i α j -1 . (3.13)
Since we wish to test the predictive ability of this relation, we split the data set into a training and a test subsets, as shown in Fig. 3.8, and we also applied a cross-validation procedure to detect an eventual bias in the choice of sets. Regression using the training set yielded c 0 = 0.25 ± 0.02 and c 1 = 0.11 ± 0.02 for the coefficients of Eq. (3.13). Its ability to predict the k ij factors from SAPT is illustrated in Fig. 3.9a, where dimers are colour-coded into four groups, depending on the charge of the monomers: cation-anion, cation-neutral, anion-neutral and neutral-neutral. The agreement is overall good, with a standard deviation of 0.07 for the training set and 0.10 for the test set, which are equivalent in terms of scatter. In order to check if a bias was introduced by our choices of training and test sets, we carried out a crossvalidation by refitting Eq. (3.13) for all the dimers leaving out each data point in turn, the so-called leave-one-out cross-validation [START_REF] Brereton | Applied Chemometrics for Scientists[END_REF]. The cross-validation did not reveal a significant bias arising from the choice of data sets.

The only systematic deviation apparent in Fig. 3.9a occurs for dimers of neutral, non-polar fragments, for which Eq. (3.13) predicts k pred = 1. For such dimers, SAPT calculations give a small induction contribution, so k SAPT < 1 with values roughly in the range 0.89 to 0.95. The deviations for dimers of non-polar fragments are commensurate with the overall scatter of the fit but a systematic trend is observed nonetheless. We attempted to improve on this small systematic deviation by considering that, even if fragments are non-polar, individual atoms can have significant partial charges, which give rise to induction effects. It could make sense to abandon the fragment approach when evaluating the second term in Eq. ( 3.13) and calculate it instead atom-by-atom,

N n M m r 2 mn (q 2 n α m + q 2 m α n ) α m α n , (3.14)
where N and M are the number of atoms in fragments i and j. This per-atom formula would require explicit interatomic distances r mn , which depend on mutual orientations of the fragments. For simplicity we replaced those by the distance between the centres of mass of the fragments, r ij . The resulting fit is shown in Fig. 3.9b. This per-atom approach indeed works better for dimers of non-polar fragments, but for all the other classes the fit becomes worse.

Since the original deviations for the non-polar dimers were small and because we are interested mainly in ionic systems, we considered the fragment-based scheme to be superior and it is the one we adopt. Thus, we propose a fragment-based scheme to predict the scaling factors needed to upgrade a fixed-charge force field into a polarisable one. This scheme requires molecular quantities such as charge, permanent dipole moment and polarisability, which have clear physical meaning and are simple to calculate using accessible quantum chemical methods.

Anticipating some of the results that will be discussed below, we found that introducing polarisation explicitly in the force field and scaling the LJ terms leads to slightly lower densities for most ionic liquids, when simulations are compared to experiments. This systematic deviation is on average about -2%, a magnitude similar to the scatter observed for different ionic liquids. So, maximum density deviations with respect to experiment reach -4% for ionic liquids we studied, as shown in Fig. 3.10. We think density is an essential property that affects the ability to predict many others, so we propose to improve the calculation of density by reducing the diameters of LJ sites, σ, by 1.5% across the specification of the force field. This scaling of σ by a factor of 0.985 corrects the bias in liquid densities. Our scheme to construct a polarisable force field for ionic liquids and eutectic solvents, CL&Pol, starting from the fixed-charge CL&P model, involves the following steps, which are illustrated in Fig. 3.11:

1. Adding Drude induced dipoles derived from atomic polarisabilities; 2. Scaling down the well depth of der Waals interactions (Lennard-Jones ) by the appropriate k ij factors, obtained either:

a) from SAPT calculations on dimers (computationally expensive), or b) from a general predictive scheme, which requires knowledge of atomic polarisabilities, fragment dipole moments and dimer geometries from straightforward quantum chemical methods;

3. Correcting the density by scaling down atomic diameters (Lennard-Jones σ) by 0.985.

CL&P add Drude

! α i scale LJ ! ϵ ! k ij SAPT2+/aDZ ! μ i scale LJ ! σ 0.985 CL&Pol MP2/cc-pVTZ ! E ind , E disp MO6-2X/PVTZ Fragment-based

calculation of equilibrium and transport properties

Structural and dynamic quantities are among the most important ones that are studied by MD simulations of ionic liquids, since the absence of volatility of these fluids complicates the study of energetic properties (vapor pressure is one of the keys leading to cohesive energy and chemical potential in organic liquids). In a previous study [START_REF] Padua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF], the influence of Drude induced dipoles and of scaling the LJ potential on the density, local structure, coordination numbers and diffusion coefficients of some ionic liquids was discussed. Here we extend the analysis to a much larger data set, with the following force field settings:

1. FixQ: fixed-charged CL&P force field; The results for several representative ionic liquids, including different cation and anion headgroups, and several lengths of alkyl chains, are given in Tab. 3.3, together with experimental data from the literature [START_REF] Chirico | Thermodynamic and thermophysical properties of the reference ionic liquid: 1-Hexyl-3-methylimidazolium bis[END_REF][START_REF] França | Thermophysical properties of ionic liquid dicyanamide (DCA) nanosystems[END_REF][START_REF] Freire | Thermophysical Characterization of Ionic Liquids Able To Dissolve Biomass[END_REF][START_REF] Harris | Revised and Extended Values for Self-Diffusion Coefficients of 1-Alkyl-3-methylimidazolium Tetrafluoroborates and Hexafluorophosphates: Relations between the Transport Properties[END_REF][START_REF] Neves | Systematic Study of the Thermophysical Properties of Imidazolium-Based Ionic Liquids with Cyano-Functionalized Anions[END_REF][START_REF] Nieto De Castro | Studies on the density, heat capacity, surface tension and infinite dilution diffusion with the ionic liquids[END_REF][START_REF] Noda | Pulsed-Gradient Spin-Echo 1H and 19F NMR Ionic Diffusion Coefficient, Viscosity, and Ionic Conductivity of Non-Chloroaluminate Room-Temperature Ionic Liquids[END_REF][START_REF] Sánchez | Density, Viscosity, and Surface Tension of Synthesis Grade Imidazolium, Pyridinium, and Pyrrolidinium Based Room Temperature Ionic Liquids[END_REF][START_REF] Simons | Influence of Zn 2+ and Water on the Transport Properties of a Pyrrolidinium Dicyanamide Ionic Liquid[END_REF][START_REF] Tokuda | Physicochemical Properties and Structures of Room Temperature Ionic Liquids. 1. Variation of Anionic Species[END_REF][START_REF] Tokuda | Physicochemical Properties and Structures of Room Temperature Ionic Liquids. 2. Variation of Alkyl Chain Length in Imidazolium Cation[END_REF][START_REF] Tokuda | Physicochemical Properties and Structures of Room-Temperature Ionic Liquids. 3. Variation of Cationic Structures[END_REF][START_REF] Zec | Physicochemical properties of (1-butyl-1-methylpyrrolydinium dicyanamide+γ-butyrolactone) binary mixtures[END_REF]. Some of the simulations with our final version of force field, KCSDrude, were performed both in the NpT and NVT ensembles to check if thermostats and barostats were introducing any artefacts into the calculated quantities.

The general trends in calculated properties are discussed herein. A decrease in density of ca. 1% is observed in the majority of liquids when Drude polarisation is introduced, becoming more accentuated when the are scaled. Exceptions are pirrolidinium ionic liquids, for which we do not see such density changes, and [C 2 C 1 im][BF 4 ], whose density increases slightly. Scaling σ as proposed above improves the agreement with experimental densities overall, which is within ±2% for the polarisable force field in versions CSDrude and KCSDrude. There are no significant differences in the agreement of densities when using k ij values from SAPT or from our predictive scheme.

Adding polarisation to the CL&P force field decreases the calculated viscosity and increases the simulated diffusion coefficients by one order of magnitude, leading to immensely improved predictions of dynamic quantities. Scaling the further fluidifies the systems, improving even more the prediction of transport properties. The small correction to σ aimed at improving density has negligible effects on the calculated transport properties, agree with experiments within the statistical errors.

The use of k ij scaling factors from SAPT or from our predictive scheme also lead to similar values for the calculated properties, validating the different choices and approximations made when formulating this scheme.

We made similar comparisons of density and predicted transport properties for a mixture of an ionic liquid with a molecular solvent, [C 2 C 1 im][OAc] + DMSO, with mole fraction x DMSO = 0.4. The simulations were performed using the FixQ and KCSDrude force field settings. Results and comparisons with literature data [START_REF] Araújo | Nucleic acid bases in 1-alkyl-3-methylimidazolium acetate ionic liquids: A thermophysical and ionic conductivity analysis[END_REF][START_REF] Freire | Thermophysical Characterization of Ionic Liquids Able To Dissolve Biomass[END_REF][START_REF] Oliveira | Influence of Different Inorganic Salts on the Ionicity and Thermophysical Properties of 1-Ethyl-3-methylimidazolium Acetate Ionic Liquid[END_REF][START_REF] Pereiro | Inorganic salts in purely ionic liquid media: the development of High Ionicity Ionic Liquids (HIILs)[END_REF][START_REF] Radhi | Macroscopic and Microscopic Study of 1-Ethyl-3-methyl-imidazolium Acetate-DMSO Mixtures[END_REF] are given in Tab. 3.4.

The density predicted with the polarisable force field for the pure ionic liquid is somewhat higher than the experimental value, whereas the dynamics is a bit too fast (higher diffusivities and lower viscosities than experimental). It should be noted that [C 2 C 1 im][OAc] is a highly hygroscopic, slightly protic ionic liquid, which may cause experimental errors, mainly due to uncontrolled water content. Proton transfer is expected to be too low to affect the results. In spite of these deviations on the pure ionic liquid, the polarisable force field predicts remarkably well the properties of the mixture with the molecular compound.

The CL&Pol force field presented here has comparable performance to two other polarisable force fields from the literature: the APPLE&P force field developed by [START_REF] Borodin | Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids[END_REF] for several families of cations and anions predicts densities within -2.7 % and 1.4 % and transport properties within a factor of roughly 1.5-2. The polarisable force field developed Yethiraj (Mc-Daniel et al., 2016, 2018a) for imidazolium ILs with several anions predicts densities within -3.7 % and +2.3 % and transport properties within approximately a factor of 1.5. The CL&Pol force field is more easily extendable to many families of ILs because it is fragment based and Units are: ρ/g cm -3 , D/10 -11 m 2 s -1 , η/mPa s. Units are: ρ/g cm -3 , D/10 -11 m 2 s -1 , η/mPa s.

ρ D + D - η ρ D + D - η [C 2 C 1 im][DCA] 303 K [C 4 C 1 im][NTf 2 ] 323 K Exp 1.
ρ D + D - η ρ D + D -D solv η [C 2 C 1 im][OAc] 298 K [C 2 C 1 im][OAc] -DMSO (f DMSO = 0.4) 298
the predictive scheme for scale factors is computationally simple. It is also more easily combined with compounds and materials described by force fields with functional forms based on LJ plus Coulomb potentials, which represent the vast majority of available models.

It is interesting to investigate the origin of the improved (faster) dynamics upon introduction of polarisation, since the induction terms added should be attractive (at least in a pairwise case). Analysis of radial and spatial distribution functions (RDF and SDF, respectively) between atomic sites showed some changes in the first-neighbour shells, with slightly more intense first peaks appearing at shorter distances, as shown in Fig. 3.12. This indicates stronger attraction between immediate neighbours, as expected. On the contrary, second-shell features in the radial distribution functions are less pronounced with the polarisable model, indicating that longer-range ordering is less marked, as seen in the isodensity contours in spatial distribution functions, shown in Fig. 3.13. The increased dynamics obtained with the polarisable force field is likely due to weaker long-range correlations, without disruption of first-shell structure.

Using scaled-down charges with a nonpolarisable force field leads to first-neighbour peaks that have lower intensities and are displaced to longer distances (related to the lower densities obtained). Polarisable force fields represent interactions more faithfully, when compared to the simply scaling down ionic charges used to improve dynamics in fixed-charge models, as has been pointed out by several authors (McDaniel et al., 2018b;[START_REF] Schröder | Comparing reduced partial charge models with polarizable simulations of ionic liquids[END_REF].

To sum up

We propose a transferable, polarisable force field for ionic liquids and their mixtures with molecular compounds, CL&Pol, that is a result of upgrading the widely used CL&P fixed-charged force field. In order to compensate for the addition of explicit polarisation in the form of Drude induced dipoles, the Lennard-Jones parameters are rescaled, and the scaling factor can be evaluated, either from computationally expensive SAPT quantum calculations or from a fast and general predictive scheme. Comparison with experiment showed much improved predictions of transport properties, which were poorly described by the fixed-charge model.

4

E X P L O R I N G F R E E V O L U M E O F S I -B A S E D B R A N C H E D C AT I O N I L S

Context

ILs have unique properties as solvents, which are easily tunable via side-chain functionalization or ion replacement. Here, we will show how the CL&Pol force field can be extended to new fragments and ions in the context of exploring solubilization phenomenon in new ionic liquids with branched and Si-functionalized cations.

introduction Ionic liquids are promising solvents for the solubilization of a broad range of solutes thanks to a combination of polarity, arising from the ionic nature, and hydrophobicity dut to nonpolar side-chains. Moreover, after a chemical transformation carried out in the IL [START_REF] Hallett | Room-Temperature Ionic Liquids: Solvents for Synthesis and Catalysis. 2[END_REF][START_REF] Welton | Room-Temperature Ionic Liquids. Solvents for Synthesis and Catalysis[END_REF], products may be separated or extracted either by using a non-polar solvent or by simply subjecting the mixture to reduced pressure. This gives new possibilities to carry out organic syntheses involving gaseous compounds.

Like any thermodynamic process, physical gas absorption is determined by enthalpic and entropic contributions,

∆ solv G = ∆ solv H -T ∆ solv S, (4.1)
where the former emerges from non-covalent interactions between the solute and the solvent, and the latter from the accessible free volume within the solvent, cavity shape, and from reorganization effects. ILs propose diverse intermolecular interactions including π-stacking, H-bonding, Coulomb and van der Waals forces, and nanometer-scale structuring with domains of different polarity (Canongia Lopes et al., 2006b). But their free volume is much smaller than that of molecular solvents and water. Indeed, the acceptance probability that a sphere of a certain radius placed at random in the system does not overlap any of its ions or molecules is three times higher for n C 6 F 14 than for et al., 2010;[START_REF] Gomes | Interactions of Carbon Dioxide with Liquid Fluorocarbons[END_REF]. Modification of the chemical structure of ILs, namely by side-chain functionalization, seems to be a promising way to increase free volume.

[C 4 C 1 im][PF 6 ] (Almantariotis
It has been recently shown that branching cationic side-chains with incorporated alkylsilane and siloxane groups leads to low-density and low-viscosity fluids that solubilize argon better than their linear counterparts [START_REF] Bakis | Maximizing the free volume in ionic liquids: experimental and theoretical investigation of branched cation ionic liquids[END_REF]. Substitution of one carbon by one silicon atom, and of a methylene group by an oxygen atom lower the bond rotation barriers by 10-15 kJ mol -1 providing ions with higher conformational flexibility [START_REF] Philippi | Targeted modifications in ionic liquids -from understanding to design[END_REF]. This may be related with a drop in viscosity (by 2-6 times!). In spite of the bigger atomic weight of silicon, density also decreases, and this molar volume expansion is caused by the greater bond lengths. These relatively low liquid densities suggest looser ion packing and larger free volume [START_REF] Endo | Fast solute diffusivity in ionic liquids with silyl or siloxane groups studied by the transient grating method[END_REF], which should enhance solubility of small, weakly interacting gases in these ILs.

In order to investigate the link between the new IL structures and free volume, we focus on local environments, free volume analysis and prediction of solubility trends in the new ILs with alkylsilane and siloxane side groups. We chose argon as a test solute, because this simple monoatomic gas that can explore voids in the IL without perturbing significantly the structure of the solvent.

simulation details

The systems studied belong to the imidazolium family with alkyl, alkylsilane and siloxane side-chains, combined with bis(trifluoromethanesulfonyl)imide (NTf - 2 ), dicyanamide (DCA -), tetracyanoborate (TCB -) and triflate (OTf -) anions, which the structural formulae are represented in Fig. 4.1. Argon was used as solute. The CL&Pol polarisable force field was used for modelling the cations with alkyl side-chains, and also NTf - 2 , OTf -and DCA -anions. It was extended to branched alkylsilane cationic sidechains, to the TCB -anion and to argon, according to the protocol described in Ch. 3. The fixed-charge CL&P-based force field of SiCmim + and SiOSiCmim + proposed by Castner [START_REF] Wu | Structure of ionic liquids with cationic silicon-substitutions[END_REF] was transformed into a polarisable version by adding Drude induced dipoles on each non-hydrogen atom and scaling the parameters of Lennard-Jones potential according to Eq. (3.13) to avoid double counting of induction effects.

The atomic partial charges of SiCmim + , (SiC) 2 mim + and SiOSiCmim + cations were computed using CHelpG/MP2/cc-pVTZ [START_REF] Breneman | Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis[END_REF] on previously optimized geometries in Gaussian [START_REF] Frisch | Gaussian~16 Revision B.01. Gaussian Inc[END_REF]. The coefficients in the cosine series, Eq. ( 2.10), of the Si C Si C dihedral were found by fitting the difference between the ab-initio potential energy profile of bis(trimethylsilyl)methane in the gas phase calculated at the MP2/cc-pVTZ//HF/6-31G(d) level and the MD energy profile of an isolated molecule with the target dihedral coefficients set to zero (Fig. 4.2). A detailed description of this procedure was presented in the articles on development of OPLS-AA [START_REF] Watkins | Perfluoroalkanes: Conformational Analysis and Liquid-State Properties from ab Initio and Monte Carlo Calculations[END_REF] and CL&P (Canongia [START_REF] Lopes | CL&P: A generic and systematic force field for ionic liquids modeling[END_REF] force fields. The interactions parameters of the TCB -anion were obtained from private communication with F. Philippi (Imperial College London), and OPLS-AA [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF] served as a basis for the polarizable model of argon. Atomic polarizabilities of atoms in the ILs were taken from the recent work of Schröder (Heid et al., 2018c). The polarizabilities of neutral Si and Ar atoms were taken from a widely-used database [START_REF] Schwerdtfeger | 2018 Table of static dipole polarizabilities of the neutral elements in the periodic table[END_REF], with a contribution of 1.87 Å 3 subtracted from the atomic polarizability of silicon in order to account for its positive, close to unity partial charge in imidazolium-based cations [START_REF] Bernardes | Additive polarizabilities in ionic liquids[END_REF]. The resulting values were α Si = 3.66 Å 3 and α Ar = 1.64 Å 3 . The scaling coefficients (Tab. 4.1) applied to the original, fixed-charge Lennard-Jones parameters of the interactions between fragments were evaluated using the predictive scheme (Ch. 3). Imidazolium-based cations were represented by an 1,3-dimethylimidazolium (C 1 C 1 im + ) head group and a side-chain: neopenthyl (Np) for Npmim + and Me 4 C 5 mim + ; tetramethylsilane (TMS) for SiCmim + and (SiC) 2 mim + ; and hexamethyldisiloxane ((SiC) 2 O) for SiOSiCmim + . The linear-chain 3-methyl-1-octylimidazolium was represented by 1-ethyl-3-methylimidazolium (C 2 C 1 im + ) and butane (C 4 H 10 ) fragments to be consistent with previous work. The NTf - 2 , DCA -, TCB -, OTf -anions and argon were treated as entire fragments. The LJ diameters of atoms in ions were decreased by 1.5 % in order to improve density.

Simulation boxes that contained 300 ion pairs were used for validation of the force field. Systems composed of 8 Ar atoms plus 100 ion pairs were set up to study the solvation environements around argon. 8 atoms, which allow for reasonable statistics, were placed at the vertices of a cube with an edge length of 20 Å in the middle of the box (which has L 40 Å) and tethered by a harmonic potential of 5 kcal mol -1 to prevent Ar...Ar interactions. These systems were equilibrated for 2 ns, and then 10 or 20 ns trajectories were generated for the pure

• • • C 1 C 1 im + 0.78 C 1 C 1 im + • • • NTf - 2 0.57 TMS • • • NTf - 2 0.66 Ar • • • C 2 C 1 im + 0.80 C 1 C 1 im + • • • DCA - 0.66 TMS • • • DCA - 0.49 Ar • • • NTf - 2 0.72 C 1 C 1 im + • • • TCB - 0.62 TMS • • • TCB - 0.69 Ar • • • DCA - 0.72 C 1 C 1 im + • • • TMS 0.66 TMS • • • OTf - 0.51 Ar • • • TCB - 0.77 C 1 C 1 im + • • • Np 0.68 C 4 H 10 • • • NTf - 2 0.69 Ar • • • OTf - 0.47 C 1 C 1 im + • • • (SiC) 2 O 0.68 (SiC) 2 O • • • NTf - 2 0.64 Ar • • • TMS 1.00 C 2 C 1 im + • • • NTf - 2 0.55 (SiC) 2 O • • • DCA - 0.63 Ar • • • Np 1.00 C 2 C 1 im + • • • C 4 H 10 0.78 TMS • • • TMS 1.00 Ar • • • C 4 H 10 1.00 Np • • • NTf - 2 0.68 Np • • • Np 1.00 Ar • • • Ar 1.00 Np • • • DCA - 0.60 (SiC) 2 O • • • (SiC) 2 O 1.00
ILs and the systems with Ar, respectively. The temperature-grouped Nosé-Hoover thermostat (Goloviznina et al., 2021a;[START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF] and a barostat were used to maintain temperature at 353 K and pressure at 1 bar for structural analysis. We replaced the tgNH with a Langevin thermostat when evaluating free energy of gas solvation in ILs due to stability issues. The time step used was of 1 fs. All other parameters of the polarisable force field (Drude mass, DC-DP force constant, the damping function etc.) and simulation details (cut-off radius, PPPM precision, etc.) were identical to the ones described in the previous chapter. The viscosities of ionic liquid were evaluated from 10 ns non-equilibrium MD trajectories using the periodic perturbation method [START_REF] Hess | Determining the shear viscosity of model liquids from molecular dynamics simulations[END_REF] with applied acceleration of 0.02 × 10 -5 Å fs -2 (Ch .2). The TRAVIS software [START_REF] Brehm | TRAVIS-A free analyzer for trajectories from molecular simulation[END_REF][START_REF] Brehm | TRAVIS -A Free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories[END_REF][START_REF] Gehrke | Structure and lifetimes in ionic liquids and their mixtures[END_REF] was used to compute radial distribution functions around Ar atoms in solution; to evaluate the distributions of void spheres as a function of their radii in the pure ILs; and to determine which atomic sites are exposed to cavities. Voids with a radius of 1.9 Å were chosen for the distribution analysis.

The chemical potential of Ar in the ionic liquids was calculated using the free energy perturbation (FEP) method [START_REF] Chipot | Calculating Free Energy Differences Using Perturbation Theory[END_REF], which was enabled in LAMMPS [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF] via FEP package. The system was composed of one freely moving Ar atom surrounded by 100 ion pairs of

[SiCmim][NTf 2 ] or [(SiC) 2 mim][NTf 2 ]
. The solute interactions were progressively activated in a simulation box containing the IL: first, the van der Waals interaction with the solvent and then, the electrostatic interaction of the induced dipole with the IL. A coupling parameter, λ, was progressively increased from zero (argon present in the box but not experiencing solute-solvent interactions) to one (argon fully interacting with the solvent), successively for the Lennard-Jones and Coulomb potentials with δλ steps of 0.05 and 0.10 respectively. At each λ value the box was equilibrated for 2 ns, followed by a 2 ns FEP calculation of "creation" (λ + δλ) and "annihilation" (λ -δλ) of the solute molecule, to confirm that hysteresis is negligible. Soft-cores potentials [START_REF] Beutler | Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations[END_REF] were used to avoid singularities at λ = 0 and 1.

An example of a FEP curve is depicted in Fig. 4.3. The chemical potential difference between the initial and the final states was calculated through the relation The free energy of transfer of argon from an ideal gas state to a dilute solution, at constant temperature and pressure, was computed through a thermodynamic cycle (Fig. 4.4). Because Ar has no intramolecular modes, the residual chemical potential corresponds directly to the creation of Ar in an IL, µ res = ∆µ IL , with no other contributions were present. The mole fraction solubility was then calculated from

∆µ = -kT n-1 i=0 ln V exp -(U λ i+1 -U λ )/(kT ) λ i V λ i . (4.2) -0.5 0.0 0.5 1.0 0.0 0.2 0.4 0.6 0.8 1.0 [SiCmim][Ntf 2 ] ln(<Vexp(-(U λ i+1 -U λ i )/kT)>/<V>) λ LJ λ(0 → 1) LJ λ(1 → 0) D λ(0 → 1) D λ(1 → 0)
x = p * ρ IL RT exp - µ res RT , (4.3)
where ρ IL is the molar density of the IL solvent and p * is the partial pressure of the gas (here set to 1 bar).

results and discussion

The CL&Pol force field was validated through comparison of the predicted density and viscosity of the pure ILs against experimental values as shown in Tab. 4.2. On average, density deviations do not exceed ±2 % with the exception of [SiOSiCmim][DCA]. Viscosities were slightly underestimated.

To check the ability of the CL&Pol force field to predict solubility of gases in ionic liquids, we evaluated the residual chemical potential and the solubility of Ar in [SiCmim][NTf 2 ] and [(SiC) 2 mim][NTf 2 ], and compared it to the experimental data reported for a partial pressure of 1 bar and extrapolated to 353 K, as shown in Tab. 4.3. The absolute deviations in residual chemical potential were comparatively small, from -0.9 to -0.7 kcal mol -1 . The calculated Units are: ρ/g cm -3 , ρ dev /%, η/mPa s.

chemical potential mainly depended on the van der Waals energy required to grow a LJ sphere in the bulk IL, the induced dipole contribution being very small (Fig. 4.3). The solubility of Ar is predicted as expected on a semi-quantitative level, with the values overestimated by about 3 times, but a general trend of higher solubility of Ar in [(SiC)

2 mim][NTf 2 ] than in [SiCmim][NTf 2 ]
was confirmed. We considered that further calculations involving other ILs would not provide much relevant information due to the simplicity of solute-solvent interactions, the high computational cost (10 000-15 000 CPU hours per system) and the availability of experimental data for these systems. Therefore, in this work FEP calculations were used essentially for validation of the force field. Undoubtedly, MD simulations can provide useful information on the Ar-IL structure, allowing to establish a connection with argon solubility. We analysed local environment of the solute in two sets of ionic liquids: Npmim + , SiCmim + , Me 4 C 5 mim + , C 8 mim + , (SiC) 2 mim + and SiOSiCmim + with NTf - 2 to study the influence of the side-chain functionalization, and (SiC) 2 mim + with NTf - 2 , DCA -and TCB -to reveal the effect of anions. The radial distribution functions around Ar in the above-mentioned systems are included in Fig. 4.5.

In the cation series, we do not observe significant differences upon incorporation of a silicon atom, as shown in Fig. ] affect the ratio between these two peaks, increasing the probability of finding the side-chain of a cation near argon. This effect is most pronounced when the flexible siloxane group is introduced, the ratio g(CH 3 )/g(CF

3 ) reaching almost 2 in [SiOSiCmim][NTf 2 ]. Branching of the side-chain ([C 8 mim][NTf 2 ] versus [Me 4 C 5 mim][NTf 2 ])
does not influence the intensity of the CF 3 peak, but leads to a considerable change on that of CH 3 . However, this should not be misinterpreted as a more probable Ar-IL interaction in [C 8 Replacement of the anion leads to more prominent structural changes. In [(SiC) 2 mim][NTf 2 ], the probabilities of finding CH 3 (green curve) and CF 3 (orange curve) groups are 2.3 and 1.7, respectively. In [(SiC) 2 mim][TCB] and [(SiC) 2 mim][DCA], the intensities of the CH 3 peak reach 2.9, while the CN groups (orange curve) of the anions are almost not found around argon (g(CN) 1). In addition, the contribution of the CH 3 N group (red curve) is reduced dramatically in [(SiC) 2 mim][DCA], driven by an interaction of the positively-charged CH 3 im + CH 2 moiety with the anion (Fig. 4.6), instead of the non-polar argon. These notable changes in the argon environment can be related to a different polarity and flexibility of the anions, with a preferential presence of NTf - 2 over TCB -and DCA -. Although the anions provide the determining effect on solute-solvent interactions, another relevant contribution to gas solubility is due to free volume, which we aim to increase in this study through side-chain functionalization. Free volume can be investigated through analysis of the relative occurrence of spontaneous cavities as a function of their size. As shown in Fig. 4.7, cavities with a radius of 0.5-0.7 Å dominate in all ionic liquids. The contribution of voids with r = 1.9 Å, which is the size capable of embedding argon atoms, is relatively small. The probability of finding cavities with 1.9 Å decreases with increasing ion rigidity and shortening the non-polar side-chain along the cation series SiOSiCmim + > (SiC) 2 mim + > C 8 mim + > Me 4 C 5 mim + > SiCmim + > Npmim + . The same tendency is observed in the anion series, NTf - 2 > TCB -> DCA -, being complemented by the effects of charge delocalization and strengthening of the interaction with the imidazolium core.

The trends in the free volume probabilities in the cation and anion series are in accordance with the experimental solubilities of argon, shown in Fig.

r(void) / Å [Npmim][NTf 2 ] [SiCmim][NTf 2 ] [Me 4 C 5 mim][Ntf 2 ] [C 8 mim][NTf 2 ] [(SiC) 2 mim][NTf 2 ] [SiOSiCmim][Ntf2] 1.8 1.9 2.0 0.0 0.5 1.0 1.5 2.0 2.5 r(void) / Å [(SiC) 2 mim][DCA] [(SiC) 2 mim][TCB] [(SiC) 2 mim][NTf 2 ]
1.8 1.9 2.0 , and free volume analysis seems thus to be a reliable method to choose a good solvent for a small non-polar gas.

[Npmim][NTf 2 ] [SiCmim][NTf 2 ] [Me 4 C 5 mim][Ntf 2 ] [C 8 mim][NTf 2 ] [(SiC) 2 mim][NTf 2 ] [SiOSiCmim][NTf 2 ]
x / 10 -4

Cation

[(SiC) 2 mim][DCA] [(SiC) 2 mim][TCB] [(SiC) 2 mim][NTf 2 ]
Anion Additional information can be obtained from comparison of the local structure around the cavities with the microenvironment of Ar in these ILs, through RDF analysis. Voids with a radius of 1.9 Å, were chosen on the basis of the Lennard-Jones diameter of argon atoms. An example of a simulation box with the corresponding voids, represented by the blue spheres, is shown in Fig. 4.9. Cavities are found to be hydrophobic due to exposure of long alkyl, alkylsilane and alkylsiloxane side-chains of cations and trifluoromethyl groups of anions (Fig. , is due to a dilution effect as explained above. In all ILs, the contribution of the CH 3 N group (red curve) in the cavity environment is small, the intensity of the first peak being always close to 1.

More pronounced changes in the local structure around voids can be noticed when the anion is replaced. While NTf - 2 presents its hydrophobic CF 3 group to voids, competing with a cation, TCB -and, especially, DCA -move from the cavities to the charged domains in the IL. This effect is driven by a preferential interaction of cyano-groups with the imidazolium head group of the cation, forming weak N • • • H C hydrogen bonds (similarly to what is observed in Ar-IL systems in Fig. 4

.6).

From comparison of the RDFs, we see that the solvation environment of argon is very similar to that of cavities, with the differences being observed only in the ratio of peak intensities, namely among the anions. Indeed, trifluoromethyl groups of NTf - 2 are present in spontaneous cavities more than terminal methyl groups of cations. This is opposite to the contribution of these groups to the local structure around argon. The cyano-groups of TCB -can still be found in cavity environments, whereas in argon solutions their probability is comparable to that of CH 3 N groups. In a DCA-based IL, the anion is not present neither surrounding argon, nor around the cavities. Around cavities, the intensity of peaks due to cations is almost independent on the nature of the counterions, but this independence is not verified in Ar solutions. We confirm though simulations that argon enters non-polar voids of a suitable size in the ILs, without affecting significantly the internal structure of the solvent. The interaction between Ar and ILs is mainly dispersive, thus argon is a good probe for the free volume analysis of ILs.

To sum up

The CL&Pol force field was extended to new branched and Si-functionalized cation ionic liquids, which are promising solvents for gas solubilization. We proposed a new theoretical approach to estimate the affinity between a solute and a solvent, based on the local structure and the free volume analysis. This method allowed us to identify the solvent with the highest argon solubility, [SiOSiCmim][NTf 2 ], among the studied ILs, without need of experimental measurements. These results will be used for future studies of ethylene and acetylene absorption by Si-functionalized ILs.

O N T H E P O S S I B I L I T Y O F S P O N TA N E O U S E X F O L I AT I O N O F G R A P H E N E I N I L S

Context

We aim to demonstrate that the application of the CL&P is not limited to ionic systems, and that the proposed approach of the force field transformation can be applied to other fixed-charge models (based on OPLS-AA or similar). Moreover, we will show that our model is promising and accurate to study interfaces between nano-scale materials and ILs. We will derive a polarisable force field for graphene to evaluate how favourable is its exfoliation in ionic liquids. introduction Due to its remarkable properties such as chemical stability, high mechanical strength, optical transparency, electrical and thermal conductivity [START_REF] Novoselov | A roadmap for graphene[END_REF], graphene is an outstanding material in batteries and supercapacitors [START_REF] El-Kady | Graphene for batteries, supercapacitors and beyond[END_REF], in anti-corrosion primers [START_REF] Zhu | Mass production and industrial applications of graphene materials[END_REF] and in composite materials [START_REF] Papageorgiou | Mechanical properties of graphene and graphene-based nanocomposites[END_REF]. Among a variety of experimental techniques to produce graphene nanosheets, liquid phase exfoliation is known as a simplest and most economical method [START_REF] Kauling | The Worldwide Graphene Flake Production[END_REF]. It typically requires application of external mechanical, chemical or electrochemical energy to overcome the van der Waals forces between the layers and to produce single-to few-layer sheets stabilized in an appropriate solvent [START_REF] Parviz | Challenges in Liquid-Phase Exfoliation, Processing, and Assembly of Pristine Graphene[END_REF][START_REF] Wei | Liquid-phase exfoliation of graphite for mass production of pristine few-layer graphene[END_REF].

Atkin [START_REF] Elbourne | Molecular Resolution in situ Imaging of Spontaneous Graphene Exfoliation[END_REF] found that certain imidazolium ILs with short alkyl chains and surface energies matched to that of graphite could exfoliate graphene under ambient conditions via a purely physical room-temperature process, without requiring any form of external energy input. To explain it, he proposed [START_REF] Freeman | Ambient energy dispersion and long-term stabilisation of large graphene sheets from graphite using a surface energy matched ionic liquid †[END_REF] a four-step mechanism based on 1) adsorption of imidazolium rings by π-stacking with the conjugated system of a graphene; 2) oscillations of a step edge of the graphene flake; 3) intercalation of the IL between graphene sheets; 4) electrostatic repulsion between adsorbed cation layers on adjacent sheets. MD simulations of the IL-graphene interface and calculations of the free energy associated with peeling one graphene flake from a graphite stack in the IL medium can contribute to elucidate this proposed mechanism.

Recently, [START_REF] Bordes | Graphene in ionic liquids : interactions at interfaces, exfoliation, stabilization[END_REF] and Bordes et al. (2018a[START_REF] Bordes | Dispersion and Stabilization of Exfoliated Graphene in Ionic Liquids[END_REF] demonstrated the possibility of spontaneity of the exfoliation process in major classes of ILs through MD simulations. However, the results obtained should be called into question because of the high simulation temperature of 423 K, at which the effects of using different cations or anions were weak or absent. Thermal motion and not solvent nature seems to be responsible for exfoliation. In addition, the simulations were performed using a non-polarisable force field with ionic charges scaled by 0.8, which leads to a decrease in cohesive energy of the IL, artificially enhancing the exfoliation process.

The importance of explicit polarisation effects when modelling the graphite-solvent interface was discussed in recent studies [START_REF] Ho | Polarizability effects in molecular dynamics simulations of the graphene-water interface[END_REF][START_REF] Ho | Molecular dynamics simulation of the graphene-water interface: comparing water models[END_REF][START_REF] Pykal | Ion Interactions across Graphene in Electrolyte Aqueous Solutions[END_REF][START_REF] Schyman | Exploring Adsorption of Water and Ions on Carbon Surfaces Using a Polarizable Force Field[END_REF]. Thus, Jorgensen [START_REF] Schyman | Exploring Adsorption of Water and Ions on Carbon Surfaces Using a Polarizable Force Field[END_REF] underlines that non-polarisable force fields are unable to reproduce trends observed in the highest-level quantum mechanical methods for graphite interactions with charged species, pointing to the need of using polarisable models. Indeed, these models propose strengthening of cation-anion interactions to account better for π-stacking between graphene layers, and also with conjugated moieties in the ionic liquids, which should provide more reliable predictions.

The aim of this study is to investigate the possibility of self-exfoliation of the carbon-based nanomaterial and the role of the ionic liquid, using the CL&Pol polarisable force field. We chose to model 1-ethyl-3-methylimidazolium acetate and 1-ethyl-3-methylimidazolium triflate, the same ionic liquids in which Atkin [START_REF] Elbourne | Molecular Resolution in situ Imaging of Spontaneous Graphene Exfoliation[END_REF] observed this spontaneous process.

simulation details

The ionic liquids 1-ethyl-3-methylimidazolium acetate, [C 2 C 1 im][OAc], and 1-ethyl-3-methylimidazolium triflate, [C 2 C 1 im][OTf], (Fig. 5.1) were modelled using the CL&Pol force field (Ch. 3), while a polarisable version of the OPLS-AA force field was adopted for graphite, graphene [START_REF] Wen | Dihedral-angle-corrected registry-dependent interlayer potential for multilayer graphene structures[END_REF] and N-methylpirrolidone (NMP) [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF]. According to the fragment-based approach of CL&Pol, benzene (Bz) was taken as a representative unit of graphene, while NMP and the ions 1-ethyl-3-methylimidazolium, acetate, triflate were treated as fragments on their own. The k ij factors needed to scale fragment-fragment interactions were derived from symmetry-adapted perturbation theory calculations (SAPT2+/aDZ level, Table 5.1). The atomic diameters of sites of the ionic liquid were decreased by 1.5 % in order to improve density calculations. The simulation box was prepared with Vesta [START_REF] Momma | VESTA: a three-dimensional visualization system for electronic and structural analysis[END_REF] and Packmol [START_REF] Martínez | PACKMOL: A package for building initial configurations for molecular dynamics simulations[END_REF], and contained a stack of 4 graphene layers [START_REF] Trucano | Structure of graphite by neutron diffraction[END_REF] of 60 Å × 67 Å, periodic in the x and y directions, with a flake of 60 Å × 40 Å, periodic in x, placed on top of the stack, exposing zig-zag edges [START_REF] Kobayashi | Observation of zigzag and armchair edges of graphite using scanning tunneling microscopy and spectroscopy[END_REF], terminated with hydrogen atoms. A sufficient amount of solvent was added above the stack, forming a liquid layer with a thickness

C 2 C 1 im + • • • OAc - 0.27 Bz • • • OTf - 0.64 NMP • • • NMP 0.79 C 2 C 1 im + • • • OTf - 0.49 Bz • • • OAc - 0.37 Bz • • • Bz 0.95 C 2 C 1 im + • • • Bz 0.62 Bz • • • NMP 0.89
of at least 25 Å after equilibration. The systems contained 7450 carbon atoms in the graphite stack plus 450 ion pairs of the IL or 1000 NMP molecules, which corresponds to 30 000-40 000 particles in total (including DPs). The size of the box was increased with empty space up to 120 Å in the z direction. The graphene slab was relaxed in the xy plane at 343 K/423 K and 1 bar in the NpT ensemble using the classical Nosé-Hoover thermostat and barostat. The system was equilibrated for 2 ns with a 1 fs in timestep the NVT ensemble and a 2 ns production run was used to produce generate a for structure analysis. All other parameters concerning the polarisable force field (Drude mass, atomic polarisabilities, damping function, etc.) and the simulation details (cut-off radius, PPPM precision, etc.) were identical to the ones described in the previous chapters.

The reversible work necessary to peel a flake from the graphene stack (Fig. 5.2) was computed though a calculation of the potential of mean force (PMF) [START_REF] Bordes | Polycyclic aromatic hydrocarbons as model solutes for carbon nanomaterials in ionic liquids[END_REF][START_REF] Sresht | Quantitative Modeling of MoS2-Solvent Interfaces: Predicting Contact Angles and Exfoliation Performance using Molecular Dynamics[END_REF][START_REF] Van Eerden | Potential of mean force by thermodynamic integration: Molecular-dynamics simulation of decomplexation[END_REF]. A harmonic biasing potential of 125 kcal mol -1 , distributed over the carbon atoms of the edge of the flake was applied. The carbon atoms of the top layer of the stack under the edge of the flake, as well as the carbon atoms of the opposite side of the flake, were tethered with a harmonic potential of 10 kcal mol -1 to avoid sliding. The distance between the centres of mass of the edge row of the flake and the row of graphene atoms underneath was considered as the biasing collective variable. It was sampled between 3.0 Å and 15.0 Å with a step of 0.25 Å with 100 ps of equilibration and 200 ps of acquisition of each peeling separation. The PMF was calculated using umbrella sampling and the weighed-histogram analysis method (WHAM) [START_REF] Kumar | The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method[END_REF].

The arrangement of the ionic liquid at the surface of the flake was analysed by computing axial distribution functions (normal to the graphene plane of the stack) of selected atomic sites

g i (z) = ρ i (z) < ρ i > (5.1)
where ρ i is number density of a given atomic type i. Similarly to radial distribution functions, these distribution functions tend to 1 at long range in an isotropic fluid.

results and discussion

In order to estimate how favourable the exfoliation of the graphene is in

[C 2 C 1 im][OAc] and [C 2 C 1 im][OTf]
, we calculated the PMF profiles presented in Fig. 5.3. As a reference, we considered the exfoliation in vacuum, so curves below represent favourable solvents for graphene exfoliation, whereas curves above correspond to solvents in which the peeling process becomes more difficult than in vacuum. To validate this approach, the simulations in ILs were complemented with those in NMP, which easily intercalates between the graphene sheets due to its planarity, and is also known as a good solvent for graphene due to the match between its surface tension of NMP and the surface energy of the material (Bordes et al., 2018a). At 343 K (Fig. 5.3a), the PMFs corresponding to the exfoliation in the ionic liquid and in the molecular solvent coincide with the vacuum curve for separations smaller than 5 Å. The differences between the studied solvents become prominent at larger distances, where the NMP curve is found below the reference vacuum line and the IL curves above. The PMF of [C 2 C 1 im] [OTf] continues to follow the vacuum line being above it by 105 kJ mol -1 , and that of [C 2 C 1 im][OAc] is higher by 490 kJ mol -1 at 10 Å separation. Thus, the PMF calculations show that NMP, as expected, is a favourable solvent for graphene exfoliation, but the ILs are not, which contradicts the experimental results of Atkin [START_REF] Elbourne | Molecular Resolution in situ Imaging of Spontaneous Graphene Exfoliation[END_REF].

d / (Å) vacuum [C 2 C 1 im][OAc] [C 2 C 1 im][OTf]
A further temperature increase to 423 K (the temperature used by [START_REF] Bordes | Dispersion and Stabilization of Exfoliated Graphene in Ionic Liquids[END_REF]) does not improve the exfoliation process significantly (Fig. 5.3b). Thanks to thermal motion, the work required to peel a graphene sheet from the graphite surface is reduced both in vacuum, by 19 kJ mol -1 , and in the ILs, by 250 kJ mol -1 in [C 2 C 1 im][OAc] and by 185 kJ mol -1 in [C 2 C 1 im] [OTf], at 10 Å. Similarly to the results at lower temperature, the PMF of [C 2 C 1 im][OAc] is always found much above the vacuum line after 5 Å of separations. On contrary, the PMF of [C 2 C 1 im][OTf], initially found above, crosses below the vacuum curve at 7.5 Å and stays below it at larger separations. However, the lowering does not seem significant enough to consider this IL to be a favourable solvent.

In the simulations performed by [START_REF] Bordes | Dispersion and Stabilization of Exfoliated Graphene in Ionic Liquids[END_REF] (Fig. 8 of that paper), the trends observed with [C 2 C 1 im] [OTf] completely disagree with our results, underlying striking differences between the results obtained using the polarisable and the non-polarisable models. Whereas scaling of the ion net charges weakens the interactions between ions in the first coordination shell, the polarisable force field predicts their strengthening. While an increased affinity between graphite and the ionic liquid should ameliorate the exfoliation process, the enhanced Coulomb ion-ion and π-π graphene-graphene interactions appear as dominating in the polarisable simulations, yielding a less favourable exfoliation.

Structural analysis shows a strong ordering of ions above the graphite plane. In

[C 2 C 1 im][OAc] (Fig 5.4a
), the closest layer is formed by imidazolium cations, with N atoms at 3.7 Å from the surface. The positive charge of this first layer is compensated by a broader layer of acetate anions, with O atoms being found most likely at 5.0 Å. The following N peak is situated at 7.1 Å. The separation between the cation peaks is larger than the intramolecular distance between two N atoms of the same cation, indicating that the imidazolium rings are parallel to the surface in the first layer. This is confirmed in Fig. 5.5 (left), taken from a snapshot, with a majority of cations oriented parallel and a few found perpendicular to the surface, explaining the small peak at 5.7 Å in the distribution of N atoms. The second layer of cations is well-defined and also predominantly parallel to the surface of the graphite. The anions do not show a preferable orientation in the first coordination shell resulting in a broader first peak. The following anion peak is poorly defined, meaning a less ordered second layer.

Coordination numbers (CN) were normalized by the number of atoms of the given type in each ion in order to simplify the comparison between different ILs ). This means that the positive charge of the closest layer of cations is compensated by half of the ions forming the first layer of anions. This is a result of alternation of charged species caused by preferential adsorption of cations on the graphite surface, due to the π-stacking and dispersive interactions with non-polar side-chains of the ILs. The same effect is observed for the second electroneutral layer formed at 7.5 Å with 2.3 ions per 100 Å 2 of the graphite surface, where the total charge of the second layer of cations is balanced by the half of the anions forming the first and second layers. For

[C 2 C 1 im][OTf] (Fig 5.4b
), the anions form two layers with a small peak at 3.4 Å and an intense peak at 5.5 Å from the graphite surface, separated by a well-defined layer of cations, at 3.6 Å, in the first coordination shell. Here, the parallel orientation of imidazolium rings dominates, with a minor contribution of the cations perpendicular to the surface appearing as a small peak at 5.7 Å. The split of the anion peak into several contributions arises from the presence of parallel and perpendicular orientations of the C S bond vector with the respect to the surface, as shown in Fig. 5.5 (right). A higher intensity peak at 5.5 Å points to the predominance of the perpendicular orientation, driven by dispersive interactions between the non-polar trifluoromethyl group and graphite. In the anions oriented parallel, the atoms closer to the graphite surface appear as the small peak at 3.4 Å, while atoms that are found further contribute to the main peak at 5.5 Å. In the second coordination shell, the broad cation peak is split into two at 7.4 Å and 8.5 Å, and the anion peak is weakly marked, showing a significant disordering of this coordination shell, less ordered than in [C 2 C 1 im][OAc]. This explains why graphene exfoliation is easier in the IL with triflate anion. Our simulations confirm an adsorption of imidazolium rings by π-stacking to the graphene surface, which is in agreement with earlier simulation results [START_REF] Atilhan | Theoretical Study of Low Viscous Ionic Liquids at the Graphene Interface[END_REF][START_REF] Fedorov | Probing the neutral graphene-ionic liquid interface: insights from molecular dynamics simulations[END_REF][START_REF] Wang | Influence of ionic liquid film thickness on ion pair distributions and orientations at graphene and vacuum interfaces[END_REF]. Two prominent layers of cations are found above the graphene plane. Such strong ordering can suppress the oscillations of the step edge of the graphene flake, preventing the exfoliation process. Nevertheless, we cannot totally disprove the selfexfoliation observed experimentally under ambient conditions due to the simplifications in our model of the graphite stack (regular flake shape, no defects present, representation of the edges, etc.) and the absence of chemical interactions (like the acid-base equilibrium of acetate), which go beyond the limits of classical molecular dynamics.

To sum up

We have shown that the CL&Pol force field can be easily extended to nanomaterials, and can be used to study their interaction with ILs. We focused on the exfoliation of graphene in [OTf], which was shown to be unfavourable due to strong adsorption of imidazolium cations on the graphite surface, even upon temperature increase. The anion effect was well pronounced at both temperatures, peeling one graphene flake from a graphite stack in the IL with the triflate anion is easier than in the IL having an acetate anion.

[C 2 C 1 im][OAc] and [C 2 C 1 im]

E X T E N D I N G C L & P O L T O E L E C T R O LY T E S , P R O T I C I O N I C L I Q U I D S , A N D D E E P E U T E C T I C S O LV E N T S

Context

Hydrogen bond formation is an indispensable element of IL-water mixtures, protic ionic liquids and deep eutectic solvents. Unfortunately, straightforward application of the CL&Pol to such system results often in trajectory stability issues. Here, we propose a method to prevent the "polarisation catastrophe" and to provide a correct description of strong H-bonds in the framework of extending of the CL&Pol to these new systems.

the "polarisation catastrophe"

The CL&Pol force field brought significant improvements to the representation of structure, energetics and dynamics of aprotic ionic liquids. Its extension to some other classes of ILs and related systems is challenging, since strong electrostatic interactions may cause instability problems in the induced dipoles, resulting in the so-called "polarisation catastrophe" illustrated on Fig. 6.1 and well known in simulations using polarisable force fields [START_REF] Bedrov | Molecular Dynamics Simulations of Ionic Liquids and Electrolytes Using Polarizable Force Fields[END_REF][START_REF] Borodin | Polarizable Force Field Development and Molecular Dynamics Simulations of Ionic Liquids[END_REF][START_REF] Chowdhary | A Polarizable Force Field of Dipalmitoylphosphatidylcholine Based on the Classical drude Model for Molecular Dynamics Simulations of Lipids[END_REF][START_REF] Noskov | Molecular dynamics study of hydration in ethanol-water mixtures using a polarizable force field[END_REF][START_REF] Ponder | Current Status of the AMOEBA Polarizable Force Field[END_REF][START_REF] Ribeiro | Polarization effects in molecular dynamics simulations of glassformers Ca(NO 3 ) 2 • nH 2 O, n = 4, 6 and 8[END_REF][START_REF] Yu | Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field[END_REF]. Examples of problematic systems as far as this behaviour is concerned, include electrolytes containing small ions such as Li + ; protic ionic liquids (PIL) such as ethylammonium nitrate (EAN); deep-eutectic solvents (DES) that have strong hydrogen bonds between a H-bond donor and an organic salt; and glycols, which contain strong, persistent intramolecular H-bonds.

Strong H-bonds are established for example between H atoms bound to N or O (the H-bond donor) and negatively charged atoms (the acceptor). In the OPLS-AA force field [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF], as in many other all-stom force fields, the H atoms are represented by point charges only, embedded in the Lennard-Jones sites of the heavier atoms, at a distance of about 1 Å from their centre. Molecular compounds with the persistent intramolecular H-bonds involving OH groups, such as glycols, are challenging to simulate using polarisable force fields. The absence of a repulsive core around those H atoms and the strong electrostatic pull on the Drude dipoles belonging to the H-bond acceptor may lead to unstable trajectories, which the thermostat handling the Drude degrees of freedom is unable to correct.

Another source of the "polarisation catastrophe" may be an excessive correlation between nearby induced dipoles. This can be avoided by damping inter-and intramolecular dipoledipole interactions with a Thole damping function [START_REF] Noskov | Molecular dynamics study of hydration in ethanol-water mixtures using a polarizable force field[END_REF][START_REF] Taylor | The effect of Thole functions on the simulation of ionic liquids with point induced dipoles at various densities[END_REF][START_REF] Thole | Molecular polarizabilities calculated with a modified dipole interaction[END_REF] that represents smearing of the electron clouds, as we used in CL&Pol, as de- scribed Ch .2 and 3. This device has become a standard for polarizable simulations and is implemented in the major molecular dynamics packages, but we found it unsuitable to handle strong hydrogen-bonded systems. Some alternative strategies that have been employed to correct unstable trajectories are based either on preventing a DP from leaving its DC (a hard or reflecting wall) or on applying a restoring force when a significant displacement occurs. Thus, Roux [START_REF] Yu | Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field[END_REF] suggested introducing an additional anharmonic force into the systems containing highly-polarisable halide anions,

U hyp = k hyp (r -R 0 ) 4 , (6.1)
which is activated only beyond a certain DC-DP distance, R 0 , so that the linear polarisation response is conserved for small DC-DP elongations. Nevertheless, this requires a shorter time step of 0.5 fs that noticeably affects the total simulation "clock" time.

Later, in order to increase the time step to 1 fs, these authors proposed using a "hard wall" [START_REF] Chowdhary | A Polarizable Force Field of Dipalmitoylphosphatidylcholine Based on the Classical drude Model for Molecular Dynamics Simulations of Lipids[END_REF] to limit the maximum displacement of the DP from its DC: when a DP reaches a user-defined limit, typically of 0.2-0.25 Å, it is reflected back towards the parent site along the DP-DC bond by scaling the positions and the velocities of the DC-DP pair in question (this hard-wall scheme is implemented in GROMACS [START_REF] Lemkul | Implementation of extended Lagrangian dynamics in GROMACS for polarizable simulations using the classical Drude oscillator model[END_REF] and in OpenMM [START_REF] Huang | Molecular dynamics simulations using the drude polarizable force field on GPUs with OpenMM: Implementation, validation, and benchmarks[END_REF]). Even so, the DC-DP pair can form a dipole of a fixed length corresponding to the user-defined elongation limit, causing a gradual approach of the hydrogen atoms (for example, of an amino group) to a DP, and the simulation can break down. In his polarisable force field for urea, Yethiraj [START_REF] Jeong | A Transferable Polarizable Force Field for Urea Crystals and Aqueous Solutions[END_REF] used the "hard wall" but added small repulsive cores to the hydrogen atoms, and excluded all non-bonded intramolecular interactions in order to maintain the geometry of urea planar.

Other possible solutions are to introduce an electric field-dependent polarisability, to account for non-linear polarisation effects,

α D (E) =        α D , E E 0 , α D E 0 pE p + 1 - E E 0 p , E > E 0 , (6.2)
where E is the electric field, and p and E 0 are adjustable parameters [START_REF] Kunz | Development of a Nonlinear Classical Polarization Model for Liquid Water and Aqueous Solutions: COS/D[END_REF]; or alternatively to associate a small repulsive core to the Drude particle [START_REF] Geada | Insight into induced charges at metal surfaces and biointerfaces using a polarizable Lennard-Jones potential[END_REF], as is done in the CHARMM force field [START_REF] Brooks | CHARMM: A program for macromolecular energy, minimization, and dynamics calculations[END_REF].

The most promising way to prevent instabilities in strong H-bonded systems is to add a stiff distance-dependent damping potential that keeps charge-dipole interactions finite at short distances and can be used in conjunction with the Thole damping function. For this purpose, we adapted the Tang-Toennies (TT) damping function [START_REF] Tang | An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients[END_REF], which although developed originally for short-range damping of dispersion interactions, has been successfully adapted for use in polarisable simulations of high temperature molten salts [START_REF] Jemmer | Dipole and quadrupole polarization in ionic systems: Ab initio studies[END_REF][START_REF] Salanne | Polarization effects in ionic solids and melts[END_REF]Salanne et al., , 2012a;;[START_REF] Sharma | SEM-Drude Model for the Accurate and Efficient Simulation of MgCl2-KCl Mixtures in the Condensed Phase[END_REF][START_REF] Wilson | Polarization effects in ionic systems from first principles[END_REF].

The main focus of the present chapter is the implementation and validation of specific shortrange damping of strong charge-induced dipole interactions in order to ensure the stability of molecular dynamics trajectories in problematic systems, and thus allowing the extension of the CL&Pol force field to important new classes of compounds, such as PIL, DES, electrolytes. We considered the ammonium family of protic ionic liquids, represented by ethylammonium nitrate and triethylammonium bis(trifluoromethanesulfonyl)amide. Among many possible choices of deep eutectic solvents, we chose choline chloride and tetrabutylammonium bromide salts, which we combined with ethylene glycol or urea as a hydrogen bond donor. Solutions of a lithium salt in 1-ethyl-3-methylimidazolium bis(fluorosulfonyl)imide ionic liquid were studied as an example of a system containing a densely-charged metal ion.

simulation details

We represent cholinium, ammonium, Cl -, Br -, bis(trifluoromethansulfonyl)imide (Ntf - 2 or TFSI) and bis(fluorosulfonyl)imide (FSI) with the CL&Pol force field for ionic liquids (Ch. 3). For Li + , NO - 3 , urea and ethylene glycol, we derived a polarisable version of OPLS-AA [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF], following the same route as for CL&Pol. Atomic partial charges in nitrate, urea and ethylene glycol were computed, on geometries optimised with Gaussian [START_REF] Frisch | Gaussian~16 Revision B.01. Gaussian Inc[END_REF], using the CHelpG [START_REF] Breneman | Determining atom-centered monopoles from molecular electrostatic potentials. The need for high sampling density in formamide conformational analysis[END_REF] method at the MP2/cc-pVTZ(-f) level of theory to calculate electron densities. Atomic polarisabilities were taken from recent publications (Heid et al., 2018a,c;[START_REF] Yu | Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field[END_REF].

The non-bonded attractive energies of the original, fixed-charge force fields were scaled to avoid double counting of polarisation effects [START_REF] Padua | Resolving Dispersion and Induction Components for Polarisable Molecular Simulations of Ionic Liquids[END_REF]. The scaling factors were evaluated mostly using the fragment-based predictive scheme proposed in Ch 3. The predictive scheme avoids for case-by-case SAPT calculations, which are costly. Nonetheless, the scaling factor for EAN was derived here at the SAPT2+/aDZ level because we consider that this is a archetypal protic ionic liquid, formed by small ions and the NO - 3 fragment had not yet been described in CL&Pol. We thus opted for a first-principles calculation of the respective scaling factors.

According to the fragment approach, ethylene glycol (EG) was split into two methanol (MeOH) units, triethylammonium (N + 2220 ) was represented by trimethylammonium (N + 1110 ) and butane (C 4 H 10 ) sites, tetrabutylammonium (N + 4444 ) by tetramethylammonium (N + 1111 ) and C 4 H 10 sites, whereas cholinium (Ch + ), ethylammonium (N + 2000 ), 1-ethyl-3-methylimidazolium (C 2 C 1 im + ), urea (U), Ntf - 2 , FSI -and NO - 3 were treated as single fragments. The scaling coefficients for the Lennard-Jones parameters are reported in Tab. 6.1; the corresponding σ parameters were kept unchanged (unless mentioned below). The Lennard-Jones parameters of Li + were not scaled because of the small polarisability of this cation, α Li + = 0.03 Å 3 .

Li[FSI] [C 2 C 1 im][FSI] C 2 C 1 im + FSI - 0.66
MD simulations of periodic cubic boxes containing 10 000-15 000 atoms were performed using LAMMPS [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF] with the DRUDE [START_REF] Dequidt | Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator[END_REF] package enabled, in which we implemented the Tang-Toennies damping function (Appendix C). The time step was 1 fs and the pressure was kept at 1 bar, while different temperatures were chosen for different systems according to the availability of experimental data for comparison. After 2 ns equilibrations, 10 ns trajectories of liquid systems were generated in the NpT ensemble using conventional and temperature-grouped [START_REF] Son | Proper Thermal Equilibration of Simulations with Drude Polarizable Models: Temperature-Grouped Dual-Nosé-Hoover Thermostat[END_REF] Nosé-Hoover thermostat and barostat.

The crystal structure of EAN was modelled for 1 ns after relaxation of the box. The polarisable force field parameters and simulation setup were kept unchanged from the previous chapters.

The structure factor of liquid systems was calculated using the TRAVIS software [START_REF] Brehm | TRAVIS -A Free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories[END_REF][START_REF] Hollóczki | Triphilic Ionic-Liquid Mixtures: Fluorinated and Nonfluorinated Aprotic Ionic-Liquid Mixtures[END_REF], with atomic form factors or cross sections, f i , corresponding to X-ray and neutron scattering,

S(q) = N i=1 N j=1 x i x j f i (q)f j (q)H ij (q) N i=1 x i f i (q) 2 , (6.3)
H ij (q) = 4πρ 0 r max 0 r 2 (g ij (r) -1)
sin(qr) qr dr.

(6.4)

The Green-Kubo relation, Eq. ( 2.56), was used to evaluate the shear viscosity using the offdiagonal components of the pressure tensor from equilibrium trajectories, with the tail of the autocorrelation function smoothed using an exponential decay function to eliminate the noise at long times (Ch. 2). We will use here viscosity values from equilibrium MD to tune force field parameters, but prefer to compute viscosity from non-equilibrium MD for the final force field. The viscosity was evaluated from non-equilibrium MD simulations using the periodic perturbation method, Eq. ( 2.73), with a cosine-shaped acceleration of 0.02 × 10 -5 Å fs -2 . Einstein's relation, Eq. ( 2.54), was applied to the mean square displacements to calculate diffusion coefficients, and the Yeh-Hummer correction, Eq. ( 2.55), was used to account for finite-size effects in the final results.

6.3 force field development

Charge-induced dipole damping function

Interactions between induced dipoles are dampened at short range through Thole functions, Eq. (2.52), thus avoiding instabilities in molecular systems or in aprotic ionic liquids. However, in the systems studied in this chapter, the presence of strong interactions between induced dipoles and certain positive charges, namely "naked" H atoms involved in H-bonds or very small cations, posed additional issues to the stability of trajectories, as discussed above. The occurrence of these instabilities prompted us to also apply short-range damping (or smearing) to charge-induced dipole interactions. Our choice of short-range damping functions was guided by the literature on polarisable simulations of high temperature molten salts [START_REF] Jemmer | Dipole and quadrupole polarization in ionic systems: Ab initio studies[END_REF][START_REF] Salanne | Polarization effects in ionic solids and melts[END_REF]Salanne et al., , 2012a;;[START_REF] Sharma | SEM-Drude Model for the Accurate and Efficient Simulation of MgCl2-KCl Mixtures in the Condensed Phase[END_REF][START_REF] Wilson | Polarization effects in ionic systems from first principles[END_REF], where small monoatomic ions are common. The functional form of the Tang-Toennies damping functions [START_REF] Tang | An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients[END_REF], originally devised for dispersion interactions, is kept, as

T n (r) = 1 -c • e -br n k=0 (br) k k! , (6.5)
where b is a parameter determining the range below which the interactions are damped (Fig. 6.2). For dipole-dipole and dipole-quadrupole dispersive interactions n = 6 and 8, respectively, but for charge-dipole electrostatic interactions the sum goes to order n = 4. We chose to keep c = 1 so that T n (r) → 0 when r → 0. The Tang-Toennies damping function requires the adjustment of the b parameter that sets the range of the damping/smearing effect. In the original studies on molten salts, the b parameter was obtained through force-matching based on first-principle calculations [START_REF] Salanne | Polarization effects in ionic solids and melts[END_REF]. A detailed parameterisation for the more complex molecular ions of the present systems would imply very high computational costs. Therefore, we consider the b parameter as universal and empirical, and make use of the Tang-Toennies damping function as a safeguard preventing the "polarisation catastrophe", rather than as an accurate representation of very short-range electrostatic forces. In strongly H-bonded systems, the Tang-Toennies function (T ij , for charge-dipole short-range damping) should be used in conjunction with Thole function (S ij , for dipole-dipole damping in general polarisable simulations),

U elst (r ij ) = i j =i q i q j r ij + i j =i T ij (r ij ) q i q D,j r ij + i j =i S ij (r ij ) q D,i q D,j r ij . (6.6)
Ethylene glycol, the simplest compound forming both intra-and intermolecular hydrogen bonds and a component of many DES, was chosen as a test for a scan of the b parameter in the range 3.5-6.5. The damping function was applied to the interactions between the H atoms of the OH groups and all induced dipoles present in the system. The charge-dipole interactions between 1-2 and 1-3 neighbours were excluded, and the 1-4 interactions were scaled by 0.5 as per specification of the force field. Too small b values lead to excessive weakening of the intermolecular forces between nearest neighbours, bringing an unbalance to the non-bonded interactions and affecting the properties of the system, as shown in Tab. 6.2. Too large b values cause a 0.2-0.3 Å displacement of the oxygen's DPs from their DCs, towards the hydrogen atoms, and such strong dipoles freeze the system and increase its density. The value b = 4.5 was found to be optimum, reproducing well both equilibrium and transport properties of ethylene glycol. [START_REF] Bohne | Thermal, Conductivity, Density, Viscosity, and Prandtl-Numbers of Ethylene Glycol-Water Mixtures[END_REF][START_REF] Chandrasekhar | The spectra and the relative yield of solvated electrons produced by resonant photodetachment of iodide anion in ethylene glycol in the temperature range 296 T 453 K[END_REF] Units are: ρ/g cm -3 , ρ dev /%, D/10 -10 m 2 s -1 , η/mPa s. The Yeh-Hummer correction of diffusion coefficients was not applied at this stage. NH and tgNH denote traditional and the temperature-grouped Nosé-Hoover thermostats, respectively.

Modification of atomic diameters

The absence of Lennard-Jones sites on certain hydrogen atoms can lead in the worse cases to instabilities of the simulations, which we could mitigate through charge-dipole damping, but it can also result in significant structural changes upon introduction of polarisation. Doubly ionic H-bonds [START_REF] Hunt | Hydrogen bonding in ionic liquids[END_REF], formed between a cation and an anion, could undergo significant strengthening, "freezing" the system, which we observed in EAN. The H-bond formed between hydrogen atoms of the ammonium head group and oxygen atoms of nitrate, H N • • • O N , is shortened to 1.2 Å with an increased intensity of the first RDF peak up to 14.5 (Fig. 6.3). The peak position is comparable with the covalent N 1 H N distance and significantly reduced from the 1.8-1.9 Å obtained from ab-initio MD for the H-bond [START_REF] Bodo | A prototypical ionic liquid explored by ab initio molecular dynamics and Raman spectroscopy[END_REF][START_REF] Zahn | Structure and dynamics of the protic ionic liquid monomethylammonium nitrate ([CH 3 NH 3 ][NO 3 ]) from ab initio molecular dynamics simulations[END_REF][START_REF] Zentel | Properties of hydrogen bonds in the protic ionic liquid ethylammonium nitrate[END_REF]. Since the hydrogen atoms are just point charges embedded into the nitrogen atoms they are bound to, the repulsive potential between N 1 and O N are not sufficient to compensate the increased polarisation effects, therefore we opted to modify the repulsive diameter σ ON-N1 . This enhancement of repulsion forces between cation and anion shifts the 1st RDF peak to 1.9 Å and decreases its intensity, as illustrated in Fig. 6.3.

In order to determine the optimal σ O N -N 1 value, we evaluated the ability of the models with σ in the range 3.1-4.0 Å to predict properties of EAN correctly, as presented in Table 6.3. The densities computed with the models with σ ranges of 3.20-3.25 Å and 3.75-3.80 Å were the closest to experiment, but the observed local structures were completely different, this having a strong impact on the dynamics of a liquid phase. Thus, the final value of σ O N -N 1 = 3.75 Å was chosen because it improves structure and dynamics, and leads to transport properties in reasonable agreement with experiment, preventing "freezing" of the system. Units are: σ/Å, ρ/g cm -3 , ρ dev /%, D/10 -11 m 2 s -1 , η/mPa s. The Yeh-Hummer correction to diffusion coefficients was not applied at this stage. NH denotes Nosé-Hoover thermostat.

When modelling choline chloride (ChCl) based DES and other strongly H-bonded systems containing halide ions, the same issue has been observed. In choline chloride plus ethylene glycol (ChCl EG), the H O • • • Cl distance is notably shortened (to 1.7 Å) with a significant intensity increase (up to 38) in the first peak of the RDF, as presented in Fig. 6.4. As in the previous case, modification of σ O Cl from 3.37 Å to 3.70 Å provided us with a reasonable prediction of the structure (Fig. 6.4) and dynamics of the liquid (Tab. 6.4). Units are: ρ/g cm -3 , ρ dev /%, D/10 -11 m 2 s -1 , η/mPa s. To avoid this fine tuning of the DES force field, we revised the source of the CL&P chloride parameters, which had been obtained from crystal structures of the few available ILs more than 15 years ago (Canongia Lopes et al., 2004aLopes et al., , 2006a)). At the time this was thought to better represent the fully ionic nature of ILs, instead of sourcing parameters for ions from hydration properties. But for strongly H-bonded systems, chloride parameters from OPLS-AA, which were derived from aqueous solutions [START_REF] Chandrasekhar | Energy component analysis for dilute aqueous solutions of lithium(1+), sodium(1+), fluoride(1-), and chloride(1-) ions[END_REF], much improve the situation. Indeed, Cl -has a larger LJ diameter in OPLS-AA (4.41 Å instead of 3.65 Å in the original CL&P). Adopting OPLS-AA as a source of parameters (both and σ) for halides leads to a good prediction of equilibrium and transport properties of DES. The diffusion coefficient of the HBD in ChCl EG is improved remarkably when results using the OPLS-AA halide FF is compared to the refined CL&Pol (Tab. 6.4). In addition, a 1.5 % sigma correction is once again required to improve density, which is consistent with the earlier CL&Pol development for aprotic ILs. The same approach was successfully applied for ChCl U and [N 4444 ][Br] EG DES, which will be discussed later in this chapter.

In contrast to the above two cases, no modifications are required for the protic ionic liquid triethylammonium bis(trifluoromethanesulfonyl)amide, ( [N 2220 ][NTf 2 ]), where the ammonium hydrogen atom is shielded by the side-chains, making the close approach of the voluminous anion more difficult.

prediction of properties of pil, des and electrolytes

A variety of systems were chosen to validate the new polarisable force field: protic ionic liquids, ethylammonium nitrate (EAN) and triethylammonium bis( Both density and dynamic properties of EAN predicted by polarisable simulations using both the traditional (NH) and the temperature-grouped (tgNH) Nosé-Hoover thermostats are in good agreement with experimental data, as reported in Tab. 6.5. Non-equilibrium MD simulations with the periodic perturbation method coupled with the tgNH thermostat improve the prediction of viscosity of EAN by a factor of two.

The CL&Pol force field describes the properties of ChCl-based DES well, although the diffusion coefficients of ChCl U (1:2) seem underestimated. As this DES has a tendency to absorb water, the experimental diffusion coefficients correspond to DES-W systems (with up to several wt% of water), in which the values are larger than those of an anhydrous sample. Indeed, Mantle [START_REF] D'agostino | Molecular motion and ion diffusion in choline chloride based deep eutectic solvents studied by 1H pulsed field gradient NMR spectroscopy[END_REF], who measured the diffusion coefficients of ChCl U at 323 K, reports a viscosity of 119 mPa s, while the more recent studies of Swadzba-Kwasny [START_REF] Gilmore | Thermal Properties of Choline Chloride/Urea System Studied under Moisture-Free Atmosphere[END_REF] and Buchner [START_REF] Agieienko | Densities, Viscosities, and Electrical Conductivities of Pure Anhydrous Reline and Its Mixtures with Water in the Temperature Range[END_REF] indicate larger values of 170 mPa s at 323 K.

The transport properties of [N 2220 ][NTf 2 ] are relatively well reproduced by our model, although the density values are underestimated by 2 %. This density deviation upon introducing polarisation in pure ionic liquids was discussed in Ch. 3 and can be corrected by reducing the Table 6.5 -Experimental [START_REF] D'agostino | Molecular motion and ion diffusion in choline chloride based deep eutectic solvents studied by 1H pulsed field gradient NMR spectroscopy[END_REF][START_REF] Gilmore | Thermal Properties of Choline Chloride/Urea System Studied under Moisture-Free Atmosphere[END_REF][START_REF] Ibrahim | Physical properties of ethylene glycol-based deep eutectic solvents[END_REF][START_REF] Judeinstein | Proton Conducting Ionic Liquid Organization as Probed by NMR: Self-Diffusion Coefficients and Heteronuclear Correlations[END_REF][START_REF] Mariani | Intriguing transport dynamics of ethylammonium nitrate-acetonitrile binary mixtures arising from nano-inhomogeneity[END_REF][START_REF] Matsumoto | Preparation of room temperature ionic liquids based on aliphatic onium cations and asymmetric amide anions and their electrochemical properties as a lithium battery electrolyte[END_REF][START_REF] Russina | Structure of a Binary Mixture of Ethylammonium Nitrate and Methanol[END_REF] and calculated equilibrium and dynamic properties of PIL and DES. Units are: ρ/g cm -3 , ρ dev /%, D/10 -11 m 2 s -1 , η/mPa s. a the sigma correction of 0.985 was applied; b the experimental diffusion coefficients correspond to a system with 3.9 wt% of water [START_REF] Di Pietro | Perturbing structure and dynamics in hydrated reline with lithium chloride[END_REF].

diameter of LJ sites, σ, by 1.5 %. However, for an ammonium cation with longer side-chains in [N 4444 ][Br] EG, density is underestimated by -2.2 % even with the 1.5 % sigma correction. This density prediction is still regarded as satisfactory but a future study of the universality of the sigma scaling factor will likely be necessary. The Cl&Pol force field is able to maintain the crystal structure of solid EAN, as presented in Tab. 6.6. Most unit-cell parameters are in good agreement with experiment, except b and β, the drift of which leads to an underestimation of the density, which we considered as acceptable for a force field mainly aimed at liquids. A correct description of the micro-structure of strong H-bonded systems is a criterion of reliability with respect to the use of the Tang-Toennies damping function in the CL&Pol force field. Medium-strength hydrogen bonds are formed between the hydroxyl groups in ethylene glycol, and also between the acidic hydrogen of the cation and the oxygen atoms of the anion in [N 2220 ][NTf 2 ], as illustrated in Fig. 6.6. The H • • • A distances and the D H • • • A angles (where A is a H-bond acceptor atom and D the atom attached to the donor hydrogen) are in agreement with theoretical studies of similar systems [START_REF] Hunt | Hydrogen bonding in ionic liquids[END_REF]. The presence of stronger H-bonds in the PIL when compared to EG can be explained by the enhanced cationanion electrostatic attraction. The Cl&Pol force field predicts a linear hydrogen bond with an average distance of 1.9 Å and an angle of 160°in EAN (Figs. 6.3 and 6.7). This agrees with ab-initio MD studies of alkylammonium nitrates performed using the Car-Parrinello method [START_REF] Bodo | A prototypical ionic liquid explored by ab initio molecular dynamics and Raman spectroscopy[END_REF][START_REF] Zahn | Structure and dynamics of the protic ionic liquid monomethylammonium nitrate ([CH 3 NH 3 ][NO 3 ]) from ab initio molecular dynamics simulations[END_REF] and also with density functional-based tight-binding methods [START_REF] Zentel | Properties of hydrogen bonds in the protic ionic liquid ethylammonium nitrate[END_REF]. The linearity of the hydrogen bond was confirmed experimentally by X-ray diffraction on the crystal [START_REF] Henderson | Crystal structure of the ionic liquid EtNH3NO3-Insights into the thermal phase behavior of protic ionic liquids[END_REF].

In other studies, namely using the empirical potential structure refinement (EPSR) applied to X-ray and neutron scattering data on liquid EAN, several groups [START_REF] Hayes | Amphiphilicity determines nanostructure in protic ionic liquids[END_REF](Hayes et al., , 2013;;[START_REF] Russina | Structure of a Binary Mixture of Ethylammonium Nitrate and Methanol[END_REF][START_REF] Song | Structural Heterogeneity and Unique Distorted Hydrogen Bonding in Primary Ammonium Nitrate Ionic Liquids Studied by High-Energy X-ray Diffraction Experiments and MD Simulations[END_REF] proposed the formation of a bent H-bond, with a length of 2.2-2.4 Å and an angle of 109°(Fig. 6.7), where three hydrogen atoms of the ammonium head group share H-bonds with one oxygen atom of the anion. This bent H-bond is a result of the presence of a Lennard-Jones site of diameter 2.5 Å on each hydrogen atom in the force field obtained from EPSR. Repeating the same procedure with a model where hydrogen atoms are embedded into the neighbouring nitrogen atom (σ H = 0.8 Å), followed by MD simulations with a three-body term included, led to a linear, directional hydrogen bond with a distance of 1.9-2.0 Å [START_REF] Gontrani | The Interpretation of Diffraction Patterns of Two Prototypical Protic Ionic Liquids: a Challenging Task for Classical Molecular Dynamics Simulations[END_REF]. The ambiguous interpretation of X-ray and neutron scattering data of liquid system compromises the reliability of the potential refinement procedure. We consider the that the linear description of the hydrogen bond in EAN has stronger physical background and our model, being compatible with a directional H-bond, predicts equilibrium and dynamic properties of EAN reasonably well. Therefore, the present polarisable simulations contribute to the analysis of H-bond patterns in EAN, one of the most important protic ionic liquids, both in fundamental and applied research. We studied the liquid structure of EAN by computing static structure factors that are compared to X-ray and neutron scattering data in Fig. 6.8. The computed X-ray and neutron structure factors are both in good agreement with experimental results, with divergences observed for q values smaller than 1 Å -1 for hydrogen-rich samples, which could be explained by residual inelastic scattering of light H atoms in the experimental data, not fully removed by correction methods [START_REF] Hayes | Amphiphilicity determines nanostructure in protic ionic liquids[END_REF]. The main calculated peak at 1.71 Å -1 (Bragg spacing of 3.7 Å) corresponds to cation-anion separations [START_REF] Gontrani | The Interpretation of Diffraction Patterns of Two Prototypical Protic Ionic Liquids: a Challenging Task for Classical Molecular Dynamics Simulations[END_REF], and is in agreement with the experimental peaks near 1.6 Å -1 [START_REF] Greaves | Diversity Observed in the Nanostructure of Protic Ionic Liquids[END_REF][START_REF] Song | Structural Heterogeneity and Unique Distorted Hydrogen Bonding in Primary Ammonium Nitrate Ionic Liquids Studied by High-Energy X-ray Diffraction Experiments and MD Simulations[END_REF][START_REF] Umebayashi | Liquid Structure and the Ion-Ion Interactions of Ethylammonium Nitrate Ionic Liquid Studied by Large Angle X-Ray Scattering and Molecular Dynamics Simulations[END_REF]. In the d3-EAN sample, simulations produce a clear sharp peak at 0.64 Å -1 , corresponding to a Bragg spacing of 9.8 Å, evidence of longer range ordering in the liquid. This coincides with experimental neutron scattering peaks at 0.625 Å -1 [START_REF] Hayes | Amphiphilicity determines nanostructure in protic ionic liquids[END_REF] and 0.66 Å -1 [START_REF] Atkin | The Smallest Amphiphiles: Nanostructure in Protic Room-Temperature Ionic Liquids with Short Alkyl Groups[END_REF]. This feature also appears as a small pre-peak in X-ray scattering at 0.6 Å -1 [START_REF] Greaves | Diversity Observed in the Nanostructure of Protic Ionic Liquids[END_REF][START_REF] Song | Structural Heterogeneity and Unique Distorted Hydrogen Bonding in Primary Ammonium Nitrate Ionic Liquids Studied by High-Energy X-ray Diffraction Experiments and MD Simulations[END_REF][START_REF] Umebayashi | Liquid Structure and the Ion-Ion Interactions of Ethylammonium Nitrate Ionic Liquid Studied by Large Angle X-Ray Scattering and Molecular Dynamics Simulations[END_REF], whereas in our simulation it is split into three separate contributions at 0.42, 0.69 and 1.0 Å -1 . q / Å -1 exp calc Figure 6.8 -Experimental [START_REF] Hayes | Amphiphilicity determines nanostructure in protic ionic liquids[END_REF][START_REF] Song | Structural Heterogeneity and Unique Distorted Hydrogen Bonding in Primary Ammonium Nitrate Ionic Liquids Studied by High-Energy X-ray Diffraction Experiments and MD Simulations[END_REF] and calculated structure factor of liquid EAN at 298 K; from bottom to top: x-ray diffraction, neutron diffraction data of H-EAN, d3-EAN (deuterated ammonium group), d8-EAN (fully deuterated cation) with the respective vertical offset of 0.0, 2.0, 5.0 and 8.0.

DES differ from the previous PIL systems, in which there is no competition between the different species for the H-bonding sites. In ChCl EG, the dominant interactions are those of Cl - with H atoms of cholinium and of ethylene glycol (Fig. 6.6). The strength of these interactions decreases in the order H O (Ch + ) > H O (EG) H 1 (Ch + ) according to Fig. 6.9. In addition, and due to the high concentration of EG (2:1 mole ratio with respect to ChCl), intermolecular interactions between ethylene glycol molecules are important. The interactions between the cation and the HBD are not marked, and we did not include them in Fig. 6.9. These results are in agreement with a recent ab-initio MD study [START_REF] Alizadeh | Are There Magic Compositions in Deep Eutectic Solvents? Effects of Composition and Water Content in Choline Chloride/Ethylene Glycol from Ab Initio Molecular Dynamics[END_REF], even though the intensities of certain peaks do not coincide perfectly, as expected due from small system sizes and short AIMD trajectories, as well as obvious differences in the Hamiltonians. The positions of the first peaks in the H O Cl, H OG Cl, H 1 Cl and H OG O HG RDFs are 2.05 Å, 2.12 Å, 2.85 Å and 2.00 Å from the CL&Pol force field, and 2.10 Å, 2.14 Å, 2.75 Å and 1.80 Å from AIMD, which is quite close. The peak intensities are in good agreement between the two methods, 16.7, 7.1, 2.2 and 1.8 from the CL&Pol force field and 10.1, 7.6, 2.2 and 2.2 from AIMD, except the intensity for H O Cl, which is higher from our force field. In ChCl U, the interaction trends are similar to those of ChCl EG. The intensity of the first peak decreases in the order H O (Ch + ) > H U (U) H 1 (Ch + ) according to Fig. 6.9, and the urea-urea interaction is still noticeable due to the high concentration of this compound. The first peak positions and intensities are in a good agreement with neutron diffraction results [START_REF] Hammond | Liquid structure of the choline chlorideurea deep eutectic solvent (reline) from neutron diffraction and atomistic modelling[END_REF]: 2.06 Å, 2.31 Å, 2.88 Å and 2.15 Å from simulations and 2.10 Å, 2.20 Å, 3.03 Å and 1.95 Å from experiments for H O Cl, H U Cl, H 1 Cl and H U O U , respectively. The only significant discrepancy is observed for the H O Cl peak, which is found to be twice higher in the MD results. Nonetheless, following the previous discussion on the validation of the EAN force field, the empirical potential fitting of scattering results can lead to unreliable results.

In [N 4444 ][Br] EG (Fig. 6.10), the first peak of H OG Br is found at 2.19 Å, which agrees well with a 2.20 Å value from a recent MD study on [N 4444 ][Br] EtOH [START_REF] Takamuku | Solvation Structures of Tetraethylammonium Bromide and Tetrafluoroborate in Aqueous Binary Solvents with Ethanol, Trifluoroethanol, and Acetonitrile[END_REF]. This peak is shifted to longer distances, when compared to ChCl EG (Fig. 6.9), due to the larger size of the halide ion. A weak hydrogen bond with the cation head group protons, H 1 Br, appears at 2.84 Å, being similar to the H 1 Cl value. The interaction between EG molecules is characterized by a peak at 2.04 Å, which coincides with that in ChCl EG.

A lithium-based electrolyte, Li [FSI] [C 2 C 1 im] [FSI], was the next system of interest to validate the CL&Pol force field with the Tang-Toennies damping function. The density values of the mixture of the ionic liquid with the salt, given in (Tab. 6.7), are reproduced within 1 % for x [Li][FSI] = 0.2 and 0.4, when using the traditional NH thermostat, while the tgNH thermostat leads to slightly higher densities, by 1.5 %, which is acceptable. Viscosity values are overestimated, especially at higher lithium salt content. This behaviour can be caused by the use of the OPLS-AA model for lithium and by the approximations used during the fit of the ACF in the Green-Kubo relation. The classical NH thermostat shows better performance than the tgNH one due to its tendency to fluidize the system (caused by incorrect kinetic energy distribution), that may compensate imperfections of the force field. The lack of experimental data on diffusion coefficients does not allow an estimation of the accuracy of the calculation of the transport properties, but the general trends are in agreement with earlier work [START_REF] Bolimowska | Investigation of Li+ Cation Coordination and Transportation, by Molecular Modeling and NMR Studies, in a LiNTf2-Doped Ionic Liquid-Vinylene Carbonate Mixture[END_REF]. Units are: ρ/g cm -3 , ρ dev /%, D/10 -11 m 2 s -1 , η/mPa s.

The strong interaction between Li + and the O atom of the anion is confirmed by the radial distribution functions (Fig. 6.11). The main peak of the RDF is at 2.1 Å with an intensity of 22.7 and the corresponding coordination number is 4.5 O atoms in the Li + first solvation shell. This interaction dominates over the interaction of the anion with the H atoms of the imidazolium cation, of which the intensity of the first peak does not exceed 2. A similar behaviour of Li + was observed in the Li[Ntf 2 ] [C 6 C 1 im][Ntf 2 ] system with a peak of intensity 20 at 1.9 Å [START_REF] Bolimowska | Investigation of Li+ Cation Coordination and Transportation, by Molecular Modeling and NMR Studies, in a LiNTf2-Doped Ionic Liquid-Vinylene Carbonate Mixture[END_REF].

To sum up

We extended the polarisable CL&Pol force field, initially developed for general aprotic ionic liquids, to systems containing strong hydrogen bonds or small ions with high charge density, such as protic ionic liquids, deep eutectic solvents and alkali metal electrolytes. The Tang-Toennies damping function was used to smear the interactions between problematic charged sites and the Drude induced dipoles at short range, stabilizing the trajectories and preventing the "polarisation catastrophe".

D E S I G N I N G N E W D E S -B A S E D E L E C T R O LY T E S F O R L I T H I U M -I O N B AT T E R I E S

Context

We designed our CL&Pol force field in such a way that it can be easily mixed with existing Drude-based models, for example with the SWM4-NDP model of water. Water is one of the central compounds in nature and a key sustainable solvent. But it also can be an undesirable impurity in ILs and DES, which should be taken into account during simulations, to reproduce experimental compositions. Here we will combine the CL&Pol with the polarisable SWM4-NDP water model to study the properties of deep eutectic electrolytes for lithium-ion batteries, and the influence of small amounts of water on their properties. introduction Owing to their good rechargeable properties and efficiency, Li-ion batteries (LIBs) currently dominate on the global electrochemical energy storage market [START_REF] Armand | Lithium-ion batteries -Current state of the art and anticipated developments[END_REF]. However, their safety is not excellent since any external thermal, mechanical or electrical abuse may result in a temperature rise, followed by electrode degradation, solvent decomposition and gas formation, triggering explosion and leakage of toxic chemicals [START_REF] Jaumaux | Non-Flammable Liquid and Quasi-Solid Electrolytes toward Highly-Safe Alkali Metal-Based Batteries[END_REF].

DES are a good alternative to highly flammable organic solvents (Di Pietro et al., 2021a;[START_REF] Hansen | Deep Eutectic Solvents: A Review of Fundamentals and Applications[END_REF]. Their application in electrochemistry is growing rapidly due to their low flammability, wide liquid range, and economic advantages, since they are prepared from cheap and environmentally-friendly raw materials.

While binary mixtures of lithium salts (NTf - 2 , ClO - 4 , NO - 3 , etc.) with amides (urea, acetamide, etc.) are well covered by experimental and theoretical studies, being the subject of a recent review (Di Pietro et al., 2021a), choline chloride-based electrolytes with lithium additives remain underexplored. Important physicochemical trends in these systems, such as density and viscosity increase, and lowering of diffusion coefficients upon salt addition have been reported in the literature [START_REF] Barik | How Does Addition of Lithium Salt Influence the Structure and Dynamics of Choline Chloride-Based Deep Eutectic Solvents?[END_REF][START_REF] Dhingra | Effect of lithium chloride on the density and dynamic viscosity of choline chloride/urea deep eutectic solvent in the temperature range[END_REF][START_REF] Millia | Bio-inspired choline chloride-based deep eutectic solvents as electrolytes for lithium-ion batteries[END_REF], but some of these results suffer from poor control of water content.

To design, develop and optimize deep eutectic electrolytes for LIBs, it is important to obtain insights on the nanoscale structure of such systems, to understand the role of the lithium cations in the intermolecular network, as well as to obtain a reliable description of equilibrium and transport properties. A prototypical system composed of ChCl U and the lithium salt of the same anion, LiCl, is here the object of a joint study involving experimental Nuclear Magnetic Resonance (NMR) and polarisable molecular dynamics simulations, the results of which are reported in the present chapter.

Since DES are highly hydrophilic, capable to absorb up to 4 % of water from air in the case of ChCl U [START_REF] Chen | Water absorption by deep eutectic solvents[END_REF], the water content in the experimental sample was carefully determined and reproduced in the theoretical model. In addition, a fully anhydrous analogue, unachievable experimentally, was also considered to separate and quantify the influence of water.

simulation details

The simulated systems (Tab. 7.1) consisted of choline chloride and urea in a molar ratio of 1 : 2, plus 0.4 mol kg -1 of LiCl added, with and without 3.9 % of water, according to the experimental compositions. This water content agrees with the quantity absorbed by pure ChCl U exposed to humid air for 6 h and more [START_REF] Chen | Water absorption by deep eutectic solvents[END_REF]. Choline chloride plus urea was represented by the CL&Pol polarisable force field (Ch. 3 and 6) with the Lennard-Jones parameters of chloride taken from OPLS-AA [START_REF] Jorgensen | Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids[END_REF], and the polarisable, four-point, rigid SWM4-NDP model was used for water [START_REF] Lamoureux | A polarizable model of water for molecular dynamics simulations of biomolecules[END_REF]. The well depths of the Lennard-Jones potential were scaled down by the factors derived from the predictive scheme, Eq. (3.13) as follows. The values for DES were scaled as presented in Ch. 6, but the parameters for lithium were not modified because of its small polarisability. Since the SWM4-NDP water model was already polarisable, only the interactions of water with the CL&Pol DES fragments were scaled, and the scaling factors were determined from the molecular polarisability and the dipole moment of water,

k ij = 1 + c 1 µ 2 i α i -1 = 0.71. (7.1)
The atomic diameters of the DES components were decreased by 1.5 % in order to improve the density values. The Thole function, Eq. (2.52), was used to dampen short-range dipole-dipole interactions with a universal parameter a = 2.6. The interactions of atomic charges of lithium, naked-charge hydrogen atoms of water, cholinium and urea with Drude dipoles were smeared using the Tang-Toennies damping function, Eq. (6.5), with parameters b = 4.5 and c = 1.0. The simulation boxes were equilibrated for 2 ns, with a time step of 1 fs, and then 10 ns trajectories were generated for structure and dynamics analyses. The geometry of water molecules was kept rigid using a separate Nosé-Hoover thermostat, while the remaining species were handled using the temperature-grouped Nosé-Hoover thermostat at 353 K. The pressure of 1 bar was controlled by a Nosé-Hoover barostat. All other parameters of the polarisable force field (Drude mass, DC-DP force constant, atomic polarisabilities, etc.) and simulation details (cut-off radius, PPPM precision, etc.) are the same as in the previous chapters.

The TRAVIS software [START_REF] Brehm | TRAVIS-A free analyzer for trajectories from molecular simulation[END_REF][START_REF] Brehm | TRAVIS -A Free Analyzer and Visualizer for Monte Carlo and Molecular Dynamics Trajectories[END_REF] was used to compute the radial, angular and spatial distribution functions, and the rotation correlation times (time taken for a vector to rotate by 1 radian around any axis) of X H bonds. The diffusion coefficients were evaluated through Einstein's relation, Eq. ( 2.54), applied to the mean square displacements.The shear viscosity was evaluated from the pressure tensor using the Green-Kubo relation, Eq. (2.56); the non-equilibrium periodic perturbation method could not be considered due to technical issues related to the presence of a separate thermostat for rigid water.

results and discussion

MD simulations reveal that the strongest interaction of lithium is with the chloride anion, appearing as a high peak with a maximum at 2.32 Å, as shown in Fig. 7.1a (red curve). The peak intensity reaches 76, which is not surprising since similar probabilities were reported in the literature for lithium interacting with vinylene carbonate [START_REF] Bolimowska | Investigation of Li+ Cation Coordination and Transportation, by Molecular Modeling and NMR Studies, in a LiNTf2-Doped Ionic Liquid-Vinylene Carbonate Mixture[END_REF] and with glymes [START_REF] Saito | Li+ Local Structure in Hydrofluoroether DilutedLi-Glyme Solvate Ionic Liquid[END_REF][START_REF] Tsuzuki | Structures of [Li(glyme)]+ complexes and their interactions with anions in equimolar mixtures of glymes and Li[TFSA]: analysis by molecular dynamics simulations[END_REF]. Formed at a distance of 3.02 Å, the first coordination shell includes 3.0 chloride anions around each lithium cation, arranged in a planar triangular LiCl 2- 3 unit. The existence of similar clusters with the NTf 2 and DCA anions has been recently discussed [START_REF] Fujii | Unusual Li+ Ion Solvation Structure in Bis(fluorosulfonyl)amide Based Ionic Liquid[END_REF][START_REF] Huang | Solvation Structure and Dynamics of Li+ in Ternary Ionic Liquid-Lithium Salt Electrolytes[END_REF][START_REF] Nürnberg | Li Coordination of a Novel Asymmetric Anion in Ionic Liquid-in-Li Salt Electrolytes[END_REF]. The interaction of lithium with urea and cholinium protons, which are observed though Nuclear Overhauser Effect Spectroscopy (NOESY) [START_REF] Di Pietro | Perturbing structure and dynamics in hydrated reline with lithium chloride[END_REF], is indirect, through the chloride anion (Fig. 7.2a-d). Urea and cholinium are found in the second coordination shell of lithium at a characteristic distance of 5.05 Å (Fig. 7.1b). Although the intensity of this peak for the cholinium head group (green curve) is greater than for urea (grey curve), we should be aware of a twice smaller number of cholinium ions taken into account when normalizing the RDFs. Indeed, urea dominates in the second coordination shell with 6.75 molecules (up to r min = 7.32 Å), while only 5.6 cholinium cations are found around lithium (up to r min = 7.72 Å). Contrary to the NOE experiments, in which a quantitative detection of hydroxyl protons is problematic, MD shows the presence of 1.65 water molecules (up to r min = 5.92 Å, blue curve) in the environment of lithium, whereas the interaction with the H O group of cholinium (rose curve) is not marked. The coordination numbers presented above also include a direct interaction of lithium with the oxygen atoms of water, urea and cholinium, with a not prominent peak appearing at 2.2 Å. These interactions are suppressed by the interaction of chloride with the metal ion, since there are 10 anions per lithium cation present in the system (Tab. 7.1).

An extensive hydrogen bond network is formed by the DES components and water. Chloride appears as the major H-bond acceptor, 70 % of which interacts with hydrogen atoms of cholinium, urea and water, the remaining being captured by the electrostatic attraction of lithium. The anion forms medium strength hydrogen bonds with water (Fig. 7.3a, blue curve) and with hydroxyl groups of cholinium (rose curve) at typical distances of The oxygen atom of urea is a weaker H-bond acceptor (Fig. 7.4a and 7.2c), its interaction with H O (rose curve) and H w (blue curve) atoms appearing at 1.9 Å with a well-defined but low probability peak, while hydrogen atoms of urea (grey curve) and methyl groups of cholinium (green curve) show slightly higher intensities with their maxima at 2.15 Å and 2.58 Å, respectively. Water, on contrary, acts mainly as a H-bond donor and not as an acceptor (Fig. 7.4b): the O w . . . H interactions are less prominent, appearing as pre-peaks at 1.95 Å and 2.1 Å in the H w (blue curve) and H O (rose curve) RDFs, and no clear maximum is present at the H U (grey) curve. In addition, the peak of methyl hydrogen atoms (O w . . . H 1 ) at 2.75 Å (Fig. 7.4b, green curve) is shifted to longer distances when compared to the O U . . . H 1 RDF (Fig. 7.4a). A dominant donor role of water can be noticed from the intense peak at 3.28 Å of the O w . . . H w (Fig. 7.4b, blue curve) RDF, which represents the separation between two water molecules situated in the first coordination shell of the same chloride anion. The oxygen atom of the hydroxyl group of cholinium (Fig. 7.4b, rose curve; Fig. 7.2a) is almost not involved in H-bond interactions.

r O U ...X / Å H O (Ch + ) H 1 (Ch + ) H U (U) H w (W) 0 1 2 3 0 2 4 6 8 10 b g(r) r O w ...X / Å H O (Ch + ) H 1 (Ch + ) H U (U) H w (W)
General trends of DES-W structure are in agreement with recent ab initio [START_REF] Fetisov | First-Principles Molecular Dynamics Study of a Deep Eutectic Solvent: Choline Chloride/Urea and Its Mixture with Water[END_REF], classical MD [START_REF] Sapir | Restructuring a Deep Eutectic Solvent by Water: The Nanostructure of Hydrated Choline Chloride/Urea[END_REF] studies and with x-ray scattering experiments [START_REF] Triolo | Liquid structure of a choline chloride-water natural deep eutectic solvent: A molecular dynamics characterization[END_REF], but a detailed direct comparison is impossible due to the different water contents and to the presence of a third compound, lithium salt, in our system. The density of the system studied here is reproduced with the reasonable level of uncertainty, as can be seen in Tab. 7.2. As expected, the polarisable MD simulations predict well the diffusion coefficients, with relative deviations not exceeding 30 %. These positive deviations are reflected in minor viscosity underestimations. MD simulations complement the experimental results providing water and chloride diffusivities, the determination of which from NMR is problematic due to the significant broadening of the water peak and the extremely fast quadrupole relaxation of the chloride nuclei.

All the diffusion coefficients have the same order of magnitude, pointing to the existence of an extended intermolecular network, which affects the motion of individual species, but the exact values are not identical indicating the absence of long-lasting ion pairs. The diffusivity decreases in the order W > U > Cl -> Ch + > Li + , which does not correlate properly with size. Faster movement of urea with respect to cholinium, and even faster diffusion of water is expected from their molecular weight and observed in DES-W mixtures from the literature [START_REF] D'agostino | Molecular and ionic diffusion in aqueous -deep eutectic solvent mixtures: probing inter-molecular interactions using PFG NMR[END_REF][START_REF] Sapir | Restructuring a Deep Eutectic Solvent by Water: The Nanostructure of Hydrated Choline Chloride/Urea[END_REF]. Lithium is the smallest species in the mixture, and also the slowest one, which can be explained by its strong coordination by chloride, similarly to the behaviour observed in IL solutions [START_REF] Bolimowska | Investigation of Li+ Cation Coordination and Transportation, by Molecular Modeling and NMR Studies, in a LiNTf2-Doped Ionic Liquid-Vinylene Carbonate Mixture[END_REF][START_REF] Nürnberg | Li Coordination of a Novel Asymmetric Anion in Ionic Liquid-in-Li Salt Electrolytes[END_REF]. This coordination also lowers the diffusion coefficient of the anion when compared to the DES-W system [START_REF] Sapir | Restructuring a Deep Eutectic Solvent by Water: The Nanostructure of Hydrated Choline Chloride/Urea[END_REF] but the effect is not marked because of the significant excess of freely moving chloride. Units are: ρ/g cm -3 , ρ dev /%, D/10 -11 m 2 s -1 , η/mPa s.

Molecular reorientational times, obtained from the NMR, were compared (Tab. 7.3) to the calculated reorientational times of the H-X bond vectors, where X is an atom directly bonded to hydrogen. Although both sets of values are in good agreement, we should be aware of the simplified representation of atoms in the simulations, which does not include internal nuclear and electronic effects. The usage of bond vectors is a reasonable assumption since it takes into account the nearest atomic environment of a selected hydrogen. So, the computational study provides information about the rotation of O H groups of cholinium and water, which are difficult to obtain experimentally due to poor peak resolution. For monoatomic lithium and chloride ions, no vector can be defined, so this property cannot be determined from classical MD.

The hindrance of the rotation motion enhances in the order CH Ch < OH W < NH U < OH Ch < Li + (these are combined results of simulation and experimental studies). The shortest τ, rotation characteristic time, is observed for weakly interacting HC protons of the cholinium cation, the longest for strongly coordinated lithium, while active H-bond acceptors, namely urea, water and the hydroxyl group of cholinium, show intermediate values. Although both cholinium and water form hydrogen bonds of comparable strength and probability (Fig 7 .3 and D.1), water rotates faster thanks to it smaller size. The slowest rotation of lithium when compared to the other cations agrees with the results on binary mixtures of Li salt and ionic liquids [START_REF] Hayamizu | Studies on the translational and rotational motions of ionic liquids composed of N-methyl-N-propylpyrrolidinium (P13) cation and bis(trifluoromethanesulfonyl)amide and bis(fluorosulfonyl)amide anions and their binary systems including lithium salts[END_REF]. A fully anhydrous system can be prepared in silico to study the influence of water on the structure and dynamics of the ChCl U LiCl system. At first sight the local structure of both systems seems identical (Fig. D.2,D.3 versus Fig. 7.1,7.3,7.4), but a careful analysis reveals a decrease of the urea and cholinium contributions in the environment of lithium upon addition of water. In a water-free system, 6.0 cholinium cations and 7.4 urea molecules are found near the metal ion (at the same characteristic distances as in ChCl U LiCl W), due to enhanced Li. . . O electrostatic interactions. This can been seen in the Li. . . O U RDF in Fig. 7.5. On the contrary, lithium chloride association is unaffected, and the remaining anions form H-bonds with H O and H U atoms, which coordination numbers reaching 0.74 and 3.40 due to the absence of competition with water. These structural changes are not pronounced, with intermolecular H-bonding patterns being preserved upon water addition, typical for a DES with low water content (< 5 %) [START_REF] Shah | Effect of water on the thermo-physical properties of Reline: An experimental and molecular simulation based approach[END_REF]. More sophisticated effects, such as proton exchange between water and cholinium [START_REF] D'agostino | Molecular and ionic diffusion in aqueous -deep eutectic solvent mixtures: probing inter-molecular interactions using PFG NMR[END_REF][START_REF] Posada | Reline aqueous solutions behaving as liquid mixtures of H-bonded co-solvents: microphase segregation and formation of co-continuous structures as indicated by Brillouin and 1H NMR spectroscopies[END_REF] or rearrangements of the interaction network [START_REF] Shah | Effect of water on the thermo-physical properties of Reline: An experimental and molecular simulation based approach[END_REF] can occur when water is intentionally added to the system in quantities greater than the DES can absorb from air. Density decreases by a few percentage points when water is added to the DES mixture, as reported in Tab. 7.2. Moderate differences between ChCl U LiCl and ChCl U LiCl W systems can be noticed in the transport properties. The addition of 3.9 % of water increases diffusion coefficients by 2.5-3 times, the reorientation is accelerated twice, and viscosity drops down even with such a low water content. This behaviour is consistent with the previous works on pure DES and DES-water mixtures [START_REF] Agieienko | Densities, Viscosities, and Electrical Conductivities of Pure Anhydrous Reline and Its Mixtures with Water in the Temperature Range[END_REF][START_REF] Posada | Reline aqueous solutions behaving as liquid mixtures of H-bonded co-solvents: microphase segregation and formation of co-continuous structures as indicated by Brillouin and 1H NMR spectroscopies[END_REF], and is potentially interesting for applications in electrolytes, where low viscosity is a desirable property [START_REF] Borodin | Molecular Modeling of Electrolytes[END_REF].

To sum up

We have shown that the CL&Pol force field can be easily combined with existing polarisable models sharing a similar Drude-based architecture. We used this approach to reveal structural patterns and dynamic trends in a prototypical choline chloride deep eutectic electrolyte with lithium salt for LIBs, in a which strong coordination of lithium cation by chloride anions was observed. Addition of small amounts water has a favourable effect on the dynamic properties, decreasing the viscosity of the system while the local environment typical of anhydrous DES LiCl is only weakly affected.

C O N C L U S I O N S A N D P E R S P E C T I V E S conclusions

We developed a transferable, general, polarisable force field for aprotic ionic liquids and their mixtures with molecular compounds, CL&Pol, obtained though transformation of the widely used CL&P fixed-charged force field. The layout of the procedure followed is presented in Fig. 8.1. A fragment-based approach was applied, so that parameters for many molecular structures could be generated with ease. This helped to avoid extensive reparametrisation allowing to cover a wide variety of cation and anion structures. Exploring a multitude of molecular structures and functional groups is one of the most relevant and interesting features in ionic liquids, allowing to choose and to design alternative solvents and technological fluids with task-specific properties.

The upgrade of the fixed-charge CL&P force field requires scaling of the Lennard-Jones parameters to compensate for the addition of explicit polarisation in the form of Drude induced dipoles. The scaling of the LJ terms can be based on directly calculated factors, obtained from SAPT quantum calculations, which are expensive computationally (a few days on 16 processors). We developed an alternative predictive scheme for the scaling factors using only simple quantities such as dipole moment, polarisability and total charge of the interacting molecular fragments. This predictive scheme was validated through the calculation of density, ion diffusion coefficients and viscosity for 12 ionic liquids and one mixture with a molecular compound. Comparison with experiment showed much improved predictions of transport properties, which was a serious shortcoming of the fixed-charge model.

Analysis of the microscopic structure of the ionic liquids shows that the structure between first-neighbours is not disrupted, but the second-shell features, between ions of the same charge, are considerably weaker when polarisation is included explicitly. The faster dynamics obtained with the polarisable force field are thus linked to this loss of second-shell structure.

We extended the force field to systems containing strong hydrogen bonds or small ions with high charge density. These systems cause trajectory stability issues in polarisable simulations due to the strong pull of the small ions or of ionic H-bonds on the Drude induced dipoles. We found that a damping, or smearing, function acting at short range between those problematic charged sites and the Drude induced dipoles was needed in order to stabilize the trajectories, preventing the "polarisation catastrophe", well known in polarisable simulations. Inspired by work in the literature on molten salts, we applied a Tang-Toennies damping function to the interaction of Drude dipoles in those specific cases, and we implemented this damping function to the LAMMPS molecular dynamics code. The parameters of the damping function were validated through analyses of the structure and dynamics of several fundamentally different systems: protic ionic liquids, deep eutectic solvents, and a lithium salt in an ionic liquid electrolyte. Comparison with experiment showed good prediction of density, viscosity, ion diffusion coefficients and a correct description of the microscopic structure. The results of our polarisable simulations contribute to the understanding of the nature of hydrogen-bonding in ethylammonium nitrate, an important protic ionic liquid and intensely studied solvent, for which two conflicting views have been proposed in the literature.

We also compared the performance of traditional and temperature-grouped Nosé-Hoover thermostats for modelling ionic and ion-molecular systems. The use of the latter thermostat is strongly recommended due to the correct treatment of translational, intramolecular and polarisation degrees of freedom. At present, this thermostat is available in the OpenMM and LAMMPS molecular dynamics codes.

We recommend the periodic perturbation method of non-equilibrium MD as a method to calculate viscosity of ionic fluids, since it gives better predictions with a tighter confidence interval than using the Green-Kubo relation from equilibrium trajectories: fitting the tail of the auto-correlation function leads to large uncertainties even if long trajectories are used, and fitting of the "running integral" requires averaging over hundreds of independent trajectories.

The application of the current polarisable force field was illustrated by the simulation of an IL -nanomaterial interface, by the study of solvation of a simple monoatomic gas in ILs and by the design of a novel DES-based electrolyte for energy-storage devices. The comparison of the results obtained with a non-polarisable and a polarisable model showed a complete disagreement between the two approaches for the prediction of the graphene exfoliation trends in the IL medium, and pointed out to the importance of including explicit induction terms when modelling IL interfaces. The description of argon solvation in new imidazolium-based ILs with branched silyl and siloxyl groups contributed to the analysis of the free volume accessible to the solute in these advanced systems, providing valuable information in view of their potential use as solvents for gas capture and organic synthesis. Important structural features that considerably affect the dynamics were revealed in the local environment of a lithium salt in a choline chloride-based DES, a prototypical system used for the development of electrolytes for Li-ion batteries, and will be taken into account during the next stages of this investigation. Consequently, we demonstrated the ease of extendability of the CL&Pol model starting from OPLS-AA-based force fields without need of expensive parametrisation and its straightforward mixing with existing Drude-based polarisable models, such as SWM4-NDP water.

The fundamental result of the present work is that we have proposed a state-of-the-art tool to explore a considerable variety of ionic fluids in silico in the context of the rational design of sustainable solvents and reaction media. Our method allows a more accurate and detailed molecular description of liquid phases containing ions by providing access to target characteristics, such as affinity to specific solutes and materials. The effective cost of solvent development procedures can be significantly reduced thus accelerating the incorporation of these media into industrial processes.

perspectives

Extension to new systems

The procedure of CL&Pol development includes rescaling the Lennard-Jones potential to remove double counting of polarisation terms, for which we propose to derive the scaling factors from a general predictive scheme, to avoid computationally expensive first-principles calculations. The robustness of the prediction is demonstrated by the reliability of our approach when extending it to new structures, and the scheme was applied to dimers composed of new tetramethylphosphonium (P + 1111 ) and cholinium (Ch + ) cations, tetracyanoborate (TCB -) and 2,2,2-trifluoro-N-(trifluoromethylsulfonyl)acetamide (TSAC -) anions, dimethyl ether (Me 2 O) and tetramethylsilane (TMS) molecules. As can be seen from Fig. 8.2, our procedure remains consistent for fragments not used during the development, demonstrating its reliability then modelling systems not covered by this thesis. On the contrary, reducing atomic diameters by 1.5 % to correct an overall -2.5 % density underestimation has not been proved to be universal. For some particular cases, for example ammonium and phosphonium cations with several long alkyl side-chains, this correction is not sufficient leading to deviations of -3 to -6 %, as shown in Fig. 8.3. Therefore, the scaling factor needs to be reduced from 0.985 to 0.950 (from 1.5 % to 5 % adjustment) in order to keep the error within -2.5 to 3.0 %. This seems to be an issue arising from the very long side-chains rather than from the ionic or polar fragments. In order to confirm or reject this assumption, a study of ionic liquid crystals formed by [C 18 C 18 im] + cations with Br -, I -, DCA - and NTf - 2 anions will be soon performed in our group.

- Application and ease of access

The CL&Pol force field, developed during this thesis, has already been used in our group (Fig. 8.4) a) to model the electrical double layer formed by imidazolium-based ILs at a molybdenum disulfide electrode [START_REF] Gong | Effect of side chain modifications in imidazolium ionic liquids on the properties of the electrical double layer at a molybdenum disulfide electrode[END_REF]; b) to study the sodium diffusion in ionic liquidbased electrolytes for Na-ion batteries [START_REF] Massaro | Sodium diffusion in ionic liquid-based electrolytes for Na-ion batteries: the effect of polarizable force fields[END_REF]; c) now it is being used to investigate carbon dioxide capture, and transformation with high selectivity towards useful chemicals, in porous ionic liquids (suspensions of metal-organic frameworks particles in ILs); d) to access the structural and dynamic characteristics of thermotropic ionic liquid crystals; e) to examine the local environment and the energetics of mixtures of epoxy-functionalized IL and amine precursors for reticulated resins with charges motives.

Besides use within our group, the model gained attention of the international community, including groups leaded by Michele Parrinello (ETH Zurich, Switzerland) and Jeffrey Grossman (MIT, USA), with the corresponding map of the received requests in 2018-2021 given in Fig. 8.5. The first CL&Pol paper, published in JCTC (DOI: 10.1021/acs.jctc.9b00689) in 2019, Units are: τ/ps; R 2 is correlation coefficient, σ is the corrected sample standard deviation, where the degree of the polynomial is 4 and the cut-off is 12 Å. The pair _ coeff section of the input file requires the b parameter and the pre-exponential factor (set to 1.0 in the present work) to be specified for the interactions between the highly-charged atom and all Drude dipoles.

R 2 = 1 - N i=0 y i -f i (t i ) 2 N i=0 y i -y 2 , σ = N i=0 y i -f i (t i ) 2 N -1 . Table A.
pair _ coeff 1 2 coul/tt 4.5 1.0

This new pair style is available in the USER-DRUDE package of LAMMPS, version 29Oct2020 and later. A detailed description of the preparation of input files for strongly H-bonded systems and the usage of the coul/tt pair style is provided in the code repositories [START_REF] Di Pietro | Perturbing structure and dynamics in hydrated reline with lithium chloride[END_REF].

Résumé

Les liquides ioniques (LIs), des sels liquides à température ambiante, sont des milieux de solvatation prometteurs grâce à leurs propriétés uniques, telles qu'une faible inflammabilité, une volatilité négligeable, une conductivité élevée, et une bonne stabilité thermique et électrochimique. Les solvants eutectiques (DES), une classe de composés ayant des propriétés similaires, ne sont pas purement ioniques, étant des mélanges d'un sel avec une espèce moléculaire. Une grande variété de LIs et de DES peut être obtenue en changeant des ions, ou en modifiant des groupes fonctionnels ou les chaînes latérales, et il est coûteux de caractériser une telle gamme de composés par les techniques expérimentales. Ainsi, des méthodes théoriques sont essentielles pour la conception rationnelle de ces systèmes, et pour obtenir des informations fondamentales sur la façon dont la structure moléculaire et les interactions déterminent leurs propriétés physiques et chimiques. Leur modélisation est une tâche difficile en raison de la diversité de leurs interactions, et un bon champ de forces moléculaires est important pour fournir une description fiable des LIs et des DES.

Nous avons développé un champ de forces général, transférable et polarisable, appelé Cl&Pol, pour la simulation moléculaire de liquides ioniques protiques et aprotiques, de solvants eutectiques, et d'électrolytes. Ce modèle une amélioration majeure par rapport aux champs de forces à charge fixe existants, qui ne peuvent pas représenter avec le même niveau de réalisme physique les interactions dans les milieux ioniques et polaires, et qui échouent dans la prévision simultanée des propriétés structurelles, d'équilibre et de transport. Pour compenser l'ajout de polarisation explicite sous forme de dipôles induits de Drude, les paramètres Lennard-Jones du champ de force original ont été réduits, par un facteur d'échelle évalué soit par des calculs quantiques couteux, soit par un schéma de prédiction rapide et général que nous proposons. Des fonctions d'écran à courte portée ont été introduites dans le champ de forces, pour représenter l'écrantage des interactions à courte portée entre les charges et les dipôles induits impliquant des petits atomes fortement chargés, empêchant ainsi des instabilités dans les trajectoires. Le nouveau champ de forces donne des trajectoires stables et de bien meilleures prédictions des propriétés de transport. L'approche par fragments est une adoptée dans ce modèle : il ne s'agit pas seulement d'un modèle spécifique pour quelques composés, mais il peut être facilement élargi et combiné avec des champs de forces polarisables et non polarisables pourvu qu'ils partagent des formes fonctionnelles compatibles. Le modèle CL&Pol a été utilisé pour étudier la solvatation de colorants et de gaz dans les solvants mentionnés ci-dessus, pour décrire les interfaces des LIs avec des nanomatériaux et pour concevoir des électrolytes pour des dispositifs de stockage d'énergie.

Mots clés: polarisation, champ de forces, dynamique moléculaire, liquides ioniques, solvants eutectiques.

Abstract

Ionic liquids (ILs), salts that are liquid at room temperature, are promising solvation media due to their unique properties, such as low flammability, negligible vapour pressure, high conductivity, and good thermal and electrochemical stability. Deep eutectic solvents (DES), a related class of compounds with similar properties, are not purely ionic, being mixtures of a salt with a molecular compound. A huge variety of ILs and DES can be obtained via ion replacement, changes in functional groups or side chain modification and it is costly to characterize such a wide range of compounds by experimental techniques. So, theoretical methods are key in their rational design and to obtain fundamental information on how molecular structure and interactions determine their physical and chemical properties. Modelling of these systems is a challenging task due to the diversity of their interactions, and a good molecular force field is important to provide a reliable description of ILs and DES.

We developed a general, transferable, polarisable force field for molecular simulation of protic and aprotic ionic liquids, deep eutectic solvents, and electrolytes, named CL&Pol. This is a major upgrade from existing fixed-charge force fields, which cannot represent with the same level of physical realism the interactions in ionic and polar media, failing to predict correctly equilibrium, structural and transport properties. In order to compensate for the addition of explicit polarisation in the form of Drude induced dipoles, the Lennard-Jones parameters of the original force field are rescaled, with the scaling factor evaluated either from computationally expensive quantum calculations or from a fast and general predictive scheme that we propose. Special damping functions were introduced in the force field to represent smearing of the interactions between charges and induced dipoles at a short range involving small, highly charged atoms, thus preventing the "polarisation catastrophe". The new force field gives stable trajectories and shows much improved predictions of transport properties. The fragment approach is a major feature: this is not just a specific model for a few compounds, but a framework that can be easily extended and combined with polarisable and non-polarisable force fields sharing compatible functional forms. The CL&Pol model was used to study solvation of dyes and gases in the above-mentioned solvents, to describe the interfaces of ILs with nanomaterials, and to design electrolytes for energy-storage devices. Key words: polarisation, force field, molecular dynamics, ionic liquids, eutectic solvents.

Figure 1 . 1 -

 11 Figure 1.1 -Structural formulae of representative aprotic ionic liquid families.

Figure 1 . 2 -

 12 Figure 1.2 -Structural formulae of representative protic ionic liquid families.

Figure 1 . 3 -

 13 Figure 1.3 -Structural formulae of representative deep eutectic solvent components.

Figure 1

 1 Figure 1.4 -Multiscale modelling of ionic liquids.

Figure 1

 1 Figure 1.5 -Representation of 1-ethyl-3-methylimidazolium cation in a) quantum mechanics; b) allatomic modelling; c) united-atom modelling; d) coarse-grained simulations.

Figure 1

 1 Figure 1.6 -Schematic view of the 1-ethyl-3-methylimidazolium cation with explicit polarisation effects represented by "charge-on-spring", which describes displacement of an electronic cloud from an atomic centre.

Figure 2 .

 2 Figure 2.1 -Schematic representation of polarisable force field types based on fluctuating charges (a), point induced dipoles (b), and Drude induced dipoles (c).

Figure 2 . 2 -

 22 Figure 2.2 -A Drude induced dipole composed of a positively-charged Drude core and a negativelycharged Drude particle, connected by a harmonic spring. The mass and the charge of the Drude particle are subtracted from those of the corresponding Drude core.

Figure 2 . 3 -

 23 Figure 2.3 -Calculated temperatures T real , T COM and T Drude as a function of the simulation time of a system with set temperatures T = 298 K and T D = 1 K controlled by the tgNH thermostat.

Figure 2 .

 2 Figure 2.4 -Left: Thole damping function with a ij = 2.6 and α ij = 1.0 Å 3 . Right: short-range damping of Drude dipole-dipole interactions using a Thole function.

Figure 2 . 5 -

 25 Figure 2.5 -Potential energy (right) and density (left) of [C 2 C 1 im][DCA] at 298 K and 1 bar as a function of simulation time. The red line represents a linear fit in the range of 0.75-2.0 ns, independent of time, characteristic of a well-equilibrated system.

Figure 2 . 6 -

 26 Figure 2.6 -The mean-squared displacement of a chloride anion in ChCl U at 353 K as a function of simulation time and the corresponding linear fit in the range 2.5-9.0 ns.

Figure 2 . 8 -

 28 Figure 2.8 -Fit of the ACF tail of the p xy pressure component of [C 2 C 1 im][DCA] at 303 K in the range 0.8-36.9 ps using the exponential functions.

Figure 2 .

 2 Figure 2.10 -R 2 of the fit of the ACF tail of [C 2 C 1 im][DCA] with f 2 (t) as a function of the switch time.The maximum of R 2 appears at 0.8 ps.

Figure 2 .

 2 Figure 2.11 -Left: running integrals of the Green-Kubo relation applied to the off-diagonal pressure components of [C 2 C 1 im][DCA] at 303 K, and corresponding fits with cut-off time of 1300 ps. Right: viscosity values averaged over the fitted p xy , p yz and p zx components with standard deviations as a function of a cut-off time. The recommended cut-off time at uncertainty equal 40 % of viscosity value (Zhang et al., 2015) is indicated.

Figure 2 .

 2 Figure 2.12 -Running integrals of the Green-Kubo relation applied to the off-diagonal pressure components of [C 2 C 1 im][BF 4 ] at 323 K and corresponding fits with a cut-off time of 1300 ps evaluated from a single trajectory (left) and averaged over five independent trajectories (right).

Figure 2 .

 2 Figure 2.13 -Velocity profile of ChCl EG generated by a cosine-shaped acceleration applied along the x-axis, with periodicity along the z-axis.

Figure 2 .

 2 Figure 2.14 -Viscosities of ethylene glycol (EG, left) and ethylammonium nitrate (EAN, right) at different acceleration amplitude. The red line represents a linear fit of the data with an extrapolation to zero acceleration; the quality of the fit is evaluated with the correlation coefficient R 2 , Eq. (2.63).

Figure 3 .

 3 Figure 3.1 -Structural formulae of cations and anions of ILs.

Figure 3 . 2 -

 32 Figure 3.2 -Example of application of the fragment approach. In[C 4 C 1 im][TsO], 1-butyl-3-methylimi-dazolium is represented by 1-ethyl-3-methylimidazolium and butane sites, while tosylate is split into toluene and mesylate fragments.

  are illustrated in Fig. 3.3 with the corresponding energy values reported in Tab. 3.1. The full set of SAPT calculations is provided in Appendix B.

Figure 3 . 3 -

 33 Figure 3.3 -Energy decomposition of the interaction potential of BF - 4 with C 2 C 1 im + and with C 4 H 10 obtained from SAPT calculations. Lines correspond to calculations at the sSAPT0/jaDZ level and symbols to more accurate calculations at the equilibrium distance at the SAPT2+/aDZ level.

Figure 3

 3 Figure 3.4 -Procedure to scale the well-depth of LJ interactions in the CL&Pol force field.

Figure 3 .

 3 Figure 3.5 -Optimised geometry of the isolated C 2 C 1 im + • • • DCA -dimer showing that the anion is positioned above the plane of the imidazolium ring with no close interactions with the H atoms of the ring. The closest C H • • • N approach are observed with H 1 atoms the alkyl side groups.

Figure 3

 3 Figure 3.6 -Radial distribution functions between hydrogen atoms of the cation and terminal nitrogen atoms of the anion of [C 2 C 1 im][DCA], showing closer distances and more intense peaks involving the H CR and H CW atoms, mainly the former. Results obtained from a 10 ns MD trajectory with the SDrude force field.

Figure 3 . 8 -

 38 Figure 3.8 -Fragments constituting the dimers in the training and test sets. Although some fragments are part of both sets, the dimers considered were different.

Figure 3

 3 Figure 3.9 -Predicted k ij factors compared to values calculated from SAPT2+. Training set (empty points) and test set (filled points) obtained from fragment-based (a) and atom-based (b) approaches.

Figure 3 .

 3 Figure 3.10 -Density for several ionic liquids calculated with the polarisable model, before and after correction of the σ parameter in the Lennard-Jones interactions.

Figure 3 .

 3 Figure 3.11 -Block diagram illustrating the procedure leading to the CL&Pol force field. The blocks with grey background indicate the more expensive option to determine the k ij scaling factors from SAPT calculations. Atomic polarisabilities α derived from quantum chemical calculations were obtained from the literature (Heid et al., 2018c).

2 .

 2 ScaleQ: FixQ model with scaled charges by 0.8; 3. Drude: polarisation added to FixQ model; 4. SDrude: Drude model with scaled ij LJ by k SAPT 2+ ij ; 5. CSDrude: SDrude model with scaled LJ σ ij by 0.985; 6. KSDrude: Drude model with scaled LJ ij by k pred ij ; 7. KCSDrude: KSDrude model with scaled LJ σ ij by 0.985.

Figure 3 .

 3 Figure 3.12 -Radial distribution function between representative atoms of the cation and anion of [C 2 C 1 im][DCA] obtained with the FixQ, ScaleQ and KCSDrude force fields.

Figure 4

 4 Figure 4.1 -Structural formulae of cations and anions of ILs.

Figure 4 . 2 -

 42 Figure 4.2 -Torsion profile of the Si-C-Si-C dihedral in (SiC) 2 mim + cation. The cosine series coefficients are v 1 , v 2 , v 4 = 0, v 3 = 0.6442.

Figure 4 . 3 -

 43 Figure 4.3 -FEP calculation as a function of the activation parameter λ for argon in [SiCmim][NTf 2 ]. LJ corresponds to the activation of Lennard-Jones site, D to the Drude induced dipoles

Figure 4

 4 Figure 4.4 -Thermodynamic cycle used to calculate the residual chemical potential of Ar in an IL. The top arrow corresponds to activation of the solute in the ideal gas state (here zero because Ar is monoatomic) and the bottom arrow to activation of the solute in the IL, evaluated directly in a FEP calculation. The right arrow corresponds to the transfer of the solute from vacuum to the IL, which we aim to evaluate using this cycle.

  4.5. The terminal methyl groups of the side-chain of the cations compete with the trifluoromethyl groups of the NTf - 2 anion for the solute. The intensities of CH 3 (green curve) and CF 3 (orange curve) RDF peaks are comparable in [Npmim][NTf 2 ]. The replacement of a carbon atom by a silicon atom ([SiCmim][NTf 2 ], [(SiC) 2 mim][NTf 2 ]) or the presence of a longer a side-chain in [Me 4 C 5 mim][NTf 2

  mim][NTf 2 ] than in [Me 4 C 5 mim][NTf 2 ]. In fact, this reflects of a dilution effect due to the different number of terminal methyl groups present in the ILs (three in [Me 4 C 5 mim][NTf 2 ] and only one in [C 8 mim][NTf 2 ], Fig. 4.1), which are considered in the normalization of the RDF. The contribution of the N-methyl group of the cation (red curve) is less marked, especially for cations with a longer side-chain ([Me 4 C 5 mim][NTf 2 ], [(SiC) 2 mim][NTf 2 ]) or containing a siloxane group ([SiOSiCmim][NTf 2 ]).

Figure 4

 4 Figure 4.5 -Radial distribution function around Ar in the ILs. The green curve corresponds to carbon atoms (labelled C T or C Si ) of terminal methyl groups of the side-chain of the cations, the red curve to carbon atoms (C 1 ) of the N-methyl group of the cation, and the orange curve to carbon atoms of the anion (C BT in NTf - 2 , C N in TCB -and C ZA in DCA -), according to Fig. 4.1.

  Figure 4.6 -Radial distribution function around H CR atom of imidazolium ring of (SiC) 2 mim + in Ar plus [(SiC) 2 mim][NTf 2 ], [(SiC) 2 mim][TCB] or [(SiC) 2 mim][DCA]. The labels of the sidechain of the cation are omited.

  4.8. The only inconsistency is ob-

Figure 4 . 7 -

 47 Figure 4.7 -Probability of finding spherical cavities in ILs as a function of their radius.

Figure 4 . 8 -

 48 Figure 4.8 -Experimental argon solubility in the studied ionic liquids at 323 K.

Figure 4

 4 Figure 4.9 -Snapshot of a [Npmim][NTf 2 ] simulation box with spherical cavities (represented by blue spheres) of a radius of 1.9 Å, which is sufficient to embed argon atoms.

  4.10). In [Npmim][NTf 2 ], the probability of finding CF 3 groups of NTf - 2 (orange curve) is much larger than that of terminal CH 3 groups of the cation (green curve). The difference between the intensities of these two peaks diminishes with incorporation of silicon atoms ([SiCmim][NTf 2 ] and [(SiC) 2 mim][NTf 2 ]), with increasing length of the side-chain ([Me 4 C 5 mim][NTf 2 ]), and completely vanishes when a siloxane group is added ([SiOSiCmim][NTf 2 ]). The noticeable increase of the CH 3 peak intensity, observed in an IL with a linear side-chain, [C 8 mim][NTf 2 ], when compared to its branched analogue, [Me 4 C 5 mim][NTf 2 ]

Figure 4 .

 4 Figure 4.10 -Solvent atoms exposed to cavities in ILs. The green curve corresponds to hydrogen atoms (H C ) of terminal CH 3 groups of the side-chain, the red curve to hydrogen atoms (H 1 ) of CH 3 N group of the cation, and the orange curve to fluorine atoms of NTf - 2 or to terminal nitrogen atoms (N C and N ZA ) of cyano-groups in TCB -or DCA -anions, according to Fig. 4.1.

Figure 5

 5 Figure 5.1 -Structural formulae of the ILs and of N-methylpirrolidone.

Figure 5 . 2 -

 52 Figure 5.2 -Exfoliation of the graphene flake from the graphene stack in 1-ethyl-3-methylimidazolium triflate.

Figure 5 . 3 -

 53 Figure 5.3 -Potential of mean force of peeling graphene in different media at 343 K (a) and 423 K (b). Due to the strong mechanical oscillations of the graphene flake, the PMFs in [C 2 C 1 im][OAc] and [C 2 C 1 im][OTf] at 343 K do not converge at large separations, beyond 11 Å and 13 Å, respectively.

  . The CN of N and O atoms in [C 2 C 1 im][OAc] (Fig 5.4c) cross at 4.9 Å meaning that the first electroneutral layer is formed by 1.25 ions per 100 Å 2 of the graphite surface. This distance corresponds to the first minimum of the N axial distribution function and to the first maximum of the O one (Fig 5.4a

Figure 5 . 4 -

 54 Figure 5.4 -Axial distribution functions (a, b) and coordination numbers (c, d) of N atoms of the cation and of O atoms of the anion (normalized by the number of atoms of the given type) in [C 2 C 1 im][OAc] (a, c) and [C 2 C 1 im][OTf] (b, d), represented as a function of the distance from the top flake perpendicular to the graphite surface at 343 K. Dashed lines correspond to electroneutral layers formed at distances of intersection of coordination numbers from atoms of cation and anion.

Figure 5

 5 Figure 5.5 -Microstructure of 1-ethyl-3-methylimidazolium acetate (left) and 1-ethyl-3methylimidazolium triflate (right) near the surface of graphite at 343 K.

Figure 6

 6 Figure 6.1 -"Polarisation catastrophe" in ethylammonium nitrate. The Drude particle of oxygen atom of the anion leaves its Drude core and moves towards the hydrogen atom of the cation head group.

  Figure 6.2 -Left: Tang-Toennies damping function with c = 1 and b in the range 3.5-6.5. Right: electrostatic potential energy of two opposite unit charges without damping (black) and with damping by the Tang-Toennies function with different ranges determined by the b parameter.

Figure 6 . 3 -

 63 Figure 6.3 -Radial distribution function of O N H N in EAN with different values of the σ Lennard-Jones interaction parameter. The rose curve corresponds to the unchanged σ O N -N 1 = 3.20 Å inherited from the CL&P force field, with polarisation terms added; the blue curve corresponds to the modified Lennard-Jones parameter, σ O N -N 1 = 3.75 Å, which we considered for the CL&Pol model.

Figure 6

 6 Figure 6.4 -Radial distribution function of H O Cl in ChCl + EG with different values of the Lennard-Jones σ parameter of the chloride anion. The rose curve corresponds to the unchanged CL&P force field (σ Cl = 3.65 Å, σ O-Cl = 3.37 Å), the blue curve to the refined CL&P force field, (σ Cl = 3.65 Å, σ O-Cl = 3.70 Å) and the green curve to the OPLS-AA chloride force field (σ Cl = 4.41 Å, σ O-Cl = 3.70 Å).

  a D'Agostino et al. (2011); b CL&P ref corresponds to the refined CL&P force field explained in the text; c the sigma correction of 0.985 was applied.

  trifluoromethanesulfonyl)amide ([N 2220 ][NTf 2 ]); deep eutectic solvents, choline chrolide + urea (ChCl U, 1:2), choline chrolide + ethylene glycol (ChCl EG, 1:2), tetrabutylammonium bromide + ethylene glycol ([N 4444 ][Br] EG, 1:2); and a lithium-based electrolyte consisting of lithium bis(fluorosulfonyl)imide (Li[FSI]) and 1-ethyl-3-methylimidazolium bis(fluorosulfonyl)imide ([C 2 C 1 im][FSI]). The structural formulae of these compounds are depicted in Fig. 6.5. The calculations of equilibrium and dynamic properties were carried out at different temperatures (EAN at 260 K and 298 K, [N 2220 ][NTf 2 ] at 298 K and 340 K, ChCl U at 323 K and 350 K) or at different compositions (Li[FSI] [C 2 C 1 im][FSI] at mole fractions of lithium salt of 0.2 and 0.4), according to available experimental data for comparison.

Figure 6

 6 Figure 6.5 -Structural formulae of the compounds used in the validation of the force field.

Figure 6 . 6 -

 66 Figure 6.6 -Probability contours revealing hydrogen bonds in ethylene glycol, [N 2220 ][NTf 2 ] and ChCl EG. The x-axes represent the distances between the H atoms and the acceptor atoms. The y-axes represent the angles formed by the D H • • • A hydrogen bonds, where D is a donor atom attached to the hydrogen. The N 3 and H 3 atomic labels in the second plot correspond to H bonded to N of the ammonium head group, and O BT are O atoms of NTf - 2 .

Figure 6

 6 Figure 6.7 -Probability contours and corresponding geometries revealing hydrogen bonds in EAN for a model with (right) and without (left) the LJ sites on ammonium hydrogen atoms. The x-axes represent the distances between the H N and O N atoms. The y-axes represent the angles formed by the N 1 H N • • • O N hydrogen bonds, where D is a donor atom attached to the hydrogen. The presence of the LJ sites disturbs the linearity of N 1 H N • • • O N hydrogen bond.

H

  Figure 6.10 -Radial distribution function of Br H O , Br H 1 O HG H OG in [N 4444 ][Br] -EG at 353 K.

Figure 6

 6 Figure 6.11 -Radial distribution function of Li O BT , H CW O BT and H CR O BT in Li[FSI] [C 2 C 1 im][FSI] (x Li[FSI] = 0.2) at 338 K.

Figure 7

 7 Figure 7.1 -Radial distribution functions of selected atoms around lithium cation in the ChCl U LiCl W system with the respective vertical offset of 0.0 and 2.0.

Figure 7 . 2 -

 72 Figure 7.2 -Spatial distribution functions of selected atoms around hydroxyl (a, e) and ammonium (b, f) groups of cholinium, urea (c, g) and water (d) in the ChCl U LiCl W (top) and ChCl U LiCl (bottom) systems. Isodensity contours at a, e) 5.0 N (Ch + ), 2.7 C(U), 3.8 O (W), 60.0 Cl, 7.5 Li; b, f) 3.0 N (Ch + ), 3.6 C(U), 3.8 O(W), 6.5 Cl, 3.4 Li; c, g) 6.5 N (Ch + ), 3.6 C(U), 3.5 O (W), 12.0 Cl, 5.0 Li; d) 5.0 N (Ch + ), 2.7 C(U), 4.6 O (W), 80.0 Cl, 15.0 Li times the average density around the central ion.

Figure 7 . 4 -

 74 Figure 7.4 -Radial distribution functions of selected hydrogen atoms around oxygen atom of urea (a) and water (b) in the ChCl U LiCl W system with the respective vertical offset of 0.0 and 1.0.

Figure 7

 7 Figure 7.5 -Radial distribution functions of oxygen atoms of urea (O U ) around lithium cation in the ChCl U LiCl and ChCl U LiCl W systems.

Figure 8 .

 8 Figure 8.1 -Layout of the procedure to obtain the CL&Pol force field from the original CL&P fixedcharge model.

Figure 8 . 2 -

 82 Figure 8.2 -Predicted and SAPT-calculated k ij scaling factors. Empty points are those used in the original CL&Pol development (Ch. 3), whereas filled points are more recent involving new fragments.

Figure 8 . 3 -

 83 Figure 8.3 -Predicting the density of ILs before and after LJ σ correction by 0.985 and by 0.950. The ions [ival] -corresponds to isovalerate and [chex] -to cyclohexanoate anions.

Figure 8 .

 8 Figure 8.4 -Examples of the CL&Pol application: a) constant-potential simulation of a substituted imidazolium dicyanamide ionic liquid electrolyte near charged MoS 2 surfaces, with the polarisation of the material represented by the image charge method (Gong et al., 2021); b) modelling of the interface between a substituted phosphinium chloride ionic liquid and zeolitic imidazolate framework (ZIF-8); c) thermotropic liquid crystal formed by a substituted imidazolium bis(trifluoromethanesulfonyl)imide ionic liquid.

CL

  A M M P S I M P L E M E N TAT I O N O F T H E TA N G -T O E N N I E S D A M P I N G F U N C T I O N The Tang-Toennies damping function is implemented as a coul/tt pair style pair _ style hybrid/overlay lj/cut/thole/long 2.6 12.0 coul/tt 4 12.0

Figure D. 1 -Figure D. 2 -Figure D. 3 -

 123 Figure D.1 -Probability contours revealing hydrogen bonds in the ChCl U LiCl W system. The xaxis represents the distances between the H atoms and the chloride acceptor atoms. The y-axis represents the angles formed by the D H • • • A hydrogen bonds, where D is a donor atom attached to the hydrogen. The H O and H 1 atoms correspond to the H atoms of hydroxyl and methyl groups of choline, and H U and H w are the H atoms of urea and water, respectively. The colour bar associated with the probability of interaction is identical across all the plots.

Table 2

 2

	.1 -Atomic polarisabilities by Schröder (Heid et al., 2018a,c) used in our polarisable force field.
	Atom	α	Atom	α	Atom	α
	H	0.323 N cat	1.208 Cl	4.40
	B	0.578 N ani	1.698 Br	5.80
	C(sp 3 )	1.016 O	1.144 I	8.81
	C cat (sp 2 ) 1.122 F	0.625
	C ani (sp 2 ) 1.432 P	1.237
	C(sp)	1.587 S	1.553
	Units are: α/Å	3 .

  2.7. Integration of such ACFs can give unrealistic negative viscosity values obtained from one or several pressure components. -Autocorrelation function of the off-diagonal pressure components of [C 2 C 1 im][DCA] ionic liquid at 303 K represented at different scales. Noise accumulates after 150 ps in all components. The integrated viscosities are 19.5 mPa s, -7.2 mPa s and -7.3 mPa s (up to 10 ns) and 8.0 mPa s, 8.8 mPa s and 9.4 mPa s (up to time when the ACF crosses zero value) for p xy , p yz and p zx , respectively.

		0.6 0.8 1.0				p xy p xz p yz	0.015 0.010					p xy p xz p yz
	S ACF	0.2 0.4				S ACF	0.005
		-0.2 0.0					0.000
		0.001 0.01 -0.4	0.1	1	10	100 1000 10000	-0.005	0.1	1	10	100	1000 zero time	10000
				t / ps						t / ps
	Figure 2.7

  2.8, and consequently we excluded them. According to the correlation coefficient R 2 and standard deviation σ, listed in Appendix (Tab. A.1), the functions f 2-4 (t) perform better. Nonetheless, f 4 (t) overfits the raw ACF, uncertainties in its parameters reaches 20 %, while in f 2 (t) and f 3 (t) the deviations do not exceed 2 % accordingly to Tab. A.1.

		0.015						0.015
		0.010	0.8 ps	36.9 ps	p xy f 0 (t) f 1 (t)		0.010	0.8 ps	36.9 ps	p xy f 2 (t) f 3 (t) f 4 (t)
	ACF						ACF
	S	0.005					S	0.005
		0.000						0.000
		0.1	1	10	100	1000	10000	0.1	1	10	100	1000	10000
				t / ps						t / ps

Table 3 .

 3 1 -Dispersion and induction energies of dimers obtained with different SAPT levels at the distance corresponding to the potential energy minimum. The k ij factor corresponds to the fraction of dispersion in non-Coulomb attraction.

	Dimer	Method	E tot	E disp	E ind	k ij
	C 4 H 10 • • • BF -4	sSAPT0/jaDZ	-20.0 -11.2 -19.1 0.40
		SAPT2+/aDZ	-28.8 -20.2	16.6 0.51
	C

2 C 1 im + • • • BF -

4 sSAPT0/jaDZ -355.3 -28.1 -34.9 0.45 SAPT2+/aDZ -368.4 -46.3 -42.6 0.52 Units are: E/kJ mol -1 .

Table 3 .

 3 2 -Cohesive energy, Lennard-Jones and electrostatic contributions from MD trajectories using different force field settings. FixQ: fixed-charged CL&P force field; Drude: polarisation added to CL&P, no scaling; SDrude: Drude model with scaled ij LJ parameter by k ij from SAPT2+/aDZ.

Table 3 .

 3 3 -Properties of ionic liquids calculated using different force fields. FixQ: fixed-charge CL&P force field; ScaleQ: FixQ model with scaled charges by 0.8; Drude: FixQ model with polarisation added; SDrude: Drude model with scaled LJ by k SAPT2+ ij ; CSDrude: SDrude model with scaled LJ σ by 0.985; KSDrude: Drude model with scaled LJ by k pred ij ; KCSDrude: KSDrude model with scaled LJ σ by 0.985. Experimental values of density, ion diffusion coefficients and viscosity from the literature.

  C 1 im][BF 4 ] 323 K [C 4 C 1 r][NTf 2 ] 343 K

		100	14.0	15.0	13.9	1.414 6.60	5.20	20.6
	FixQ	1.102	1.43	1.54	362 ± 313	1.419 0.66	0.47	232 ± 43
	ScaleQ	1.057	12.8	14.6	7.6 ± 1.6	1.385 4.62	2.72	54 ± 28
	Drude	1.094	8.64	10.9	19 ± 5	1.414 2.74	1.90	58 ±
	SDrude	1.095	16.2	19.7	8.6 ± 1.4	1.375 8.66	7.86	19 ±
	CSDrude	1.123	16.6	19.4	7.6 ± 1.5	1.402 7.81	6.27	15 ±
	KSDrude	1.093	15.3	14.5	8.7 ± 1.4	1.362 11.0 10.6	11 ±
	KCSDrude NPT	1.122	14.8	18.2	8.4 ± 0.8	1.388 12.2 11.9	9.1 ± 1.6
	KCSDrude NVT	1.121	17.2	18.7	8.4 ± 1.1	1.391 11.8 10.5	14 ±
	[C 2 Exp 1.264 10.5	8.99	17.4	1.355 10.2 8.90	16.2
	FixQ	1.242	1.57	0.43	445 ± 228	1.360 0.29	0.33	554 ± 264
	ScaleQ	1.175	15.9	7.09	20 ± 2	1.326 1.90	1.47	122 ± 91
	Drude	1.228	5.53	4.25	88 ± 25	1.360 0.83	0.80	223 ± 182
	SDrude	1.251	13.2	9.38	24 ± 14	1.304 6.19	5.90	19 ±
	CSDrude	1.278	11.6	8.68	22 ± 5	1.328 7.36	6.65	15 ±
	KCSDrude NVT	1.284	13.3	10.9	20 ± 7	1.323 8.31	7.96	8.7 ± 0.5
	Exp	1.013	N/A N/A	20.0	1.000 8.87	11.6	16.1
	FixQ	1.017	1.07	0.95	208 ± 117	1.005 0.19	0.15	1061 ± 465
	ScaleQ	0.984	8.44	9.66	13.1 ± 1.5	0.967 3.45	4.93	71 ± 48
	Drude	1.009	3.36	4.03	35 ± 6	1.007 0.91	1.09	607 ± 382
	SDrude	0.983	5.37	8.76	15 ± 3	0.983 3.21	4.73	23.2 ± 1.3
	CSDrude	1.002	5.16	8.34	12.4 ± 1.2	1.003 5.21	6.14	36 ± 14
	KCSDrude NVT 1.006	6.17	9.82	23 ± 7	1.006 5.57	8.09	31 ±
			[C 4 C 1 im][DCA] 323 K	[C 4 C 1 im][BF 4 ] 343 K
	Exp	1.045	N/A N/A	13.5	1.170 8.00	8.20	18.7
	FixQ	1.049	2.12	1.98	122 ± 41	1.154 1.19	0.88	347 ± 82
	Drude	1.039	7.75	9.91	23 ± 5	1.142 7.30	5.79	82 ± 40
	SDrude	1.021	15.3	20.5	9.9 ± 0.7	1.134 9.29	13.5	19 ±
	CSDrude	1.044	16.3	20.2	7.86 ± 0.15	1.159 11.7 11.2	12 ±
			[C 6 C 1 im][NTf 2 ] 343 K	[C 6 C 1 im][BF 4 ] 353 K
	Exp	1.331	8.61	8.12	14.4	1.107 6.77	7.02	21.2
	FixQ	1.340	1.03	0.68	384 ± 435	1.093 0.97	0.67	215 ± 59
	Drude	1.336	2.87	1.91	42 ± 6	1.085 3.94	4.07	72 ± 40
	SDrude	1.291	8.99	10.1	11 ± 4	1.067 4.44	6.12	26 ±
			[C 2 C 1 im][NTf 2 ] 323 K
	Exp	1.493	11.2	7.54	15.6
	FixQ	1.494	0.75	0.41	101 ± 27
	Drude	1.488	3.58	2.33	68 ± 32
	SDrude	1.458	10.9	5.89	12 ± 3

[C 6 C 1 im][DCA] 333 K [C 4 C 1 pyrr][DCA] 323 K

Table 3 .

 3 4 -Properties of ethylmethylimidazolium acetate and its mixture with DMSO calculated using different force field settings. FixQ: fixed-charge CL&P force field; KCSDrude: FixQ model with added polarisation, scaled LJ by k pred ij and scaled LJ σ by 0.985. Experimental values taken from the literature.

  K

	Exp	1.107 0.95 0.80	115 -144	1.104 5.41 5.35 15.3	16.3
	FixQ	1.114 0.11 0.06 1659 ± 1291	1.115 0.45 0.25 1.23	609 ± 417
	KCSDrude 1.147 5.14 5.26	67 ± 20	1.119 4.96 4.90 23.7	14.1 ± 1.3

Table 4 .

 4 1 -Scaling factors for non-bonded attractive interactions evaluated though the predictive scheme

	Dimer	k pred Dimer	k pred Dimer	k pred
	Ar

Table 4 .

 4

	IL	ρ exp	ρ calc	ρ dev η exp	η calc
	[Npmim][NTf 2 ]	1.3604 1.3492 -0.82 16.6 8.5 ± 0.2
	[SiCmim][NTf 2 ]	1.3393 1.3314 -0.59 11.7 7.1 ± 0.1
		2580 1.2590	0.08 12.2 6.4 ± 0.1
	[SiOSiCmim][DCA] 0.9978 1.0224	2.47 10.4 12.0 ± 0.3
	[(SiC) 2 mim][NTf 2 ]	1.2520 1.2315 -1.64 16.1 7.9 ± 0.1
	[(SiC) 2 mim][DCA]	0.9804 0.9830	0.27 16.2 13.6 ± 0.5
	[(SiC) 2 mim][TCB]	0.9460 0.9616	1.65 12.7 25.5 ± 0.8
	[(SiC) 2 mim][OTf]	1.1301 1.1235 -0.58 31.6 16.7 ± 0.5

2 -Calculated and experimental densities and viscosities of pure ILs at 353 K [C 8 mim][NTf 2 ] 1.2729 1.2577 -1.20 13.4 8.6 ± 0.2 [Me 4 C 5 mim][NTf 2 ] 1.2778 1.2689 -0.69 38.0 13.5 ± 0.3 [SiOSiCmim][NTf 2 ] 1.

Table 4 .

 4 , and predicted residual chemical potential and argon mole fraction solubility in ILs.

	IL	µ res	x
	[SiCmim][NTf 2 ]	Exp 1.99 6.7
		Sim 1.23 18.0
	[(SiC) 2 mim][NTf 2 ] Exp 1.85 10.2
		Sim 1.02 33.3
	Units are: µ/kcal mol -1 , x/10 -4 .

Table 5 .

 5 1 -Scaling factors for non-bonded attractive interactions evaluated at SAPT2+/aDZ level

	Dimer	k SAPT Dimer	k SAPT Dimer	k SAPT

Table 6 .

 6 1 -Scaling coefficients for the modification of non-bonded attractive interactions in PIL and DES.

	System	Fragment i Fragment j	k
	Ethylene glycol	MeOH	MeOH	0.87
	ChCl EG/U	Ch +	Cl -	0.48
		Ch +	MeOH	0.66
		Cl -	MeOH	0.60
		MeOH	MeOH	0.87
		Ch +	U	0.56
		Cl -	U	0.48
		U	U	0.64
	EAN	N + 2000	NO -3	0.20
	[N 2220 ][NTf 2 ]	N + 1110	C 4 H 10	0.60
		N + 1110	NTf -2	0.50
		C 4 H 10	NTf -2	0.69
		C 4 H 10	C 4 H 10	1.00
	[N 4444 ][Br] EG	N + 1111	C 4 H 10	0.64
		N + 1111	Br -	0.50
		N + 1111	MeOH	0.65
		C 4 H 10	Br -	0.65
		C 4 H 10	MeOH	0.93
		Br -	MeOH	0.64

Table 6 .

 6 2 -Experimental

  and calculated properties of ethylene glycol at 298 K, for different values of the b parameter in the Tang-Toennies damping function.

		b	ρ	ρ dev	D	η eq	η noneq
	Exp		1.110		0.9	16.63
	FixQ		1.071 -3.5 2.8	5 ± 1
	Drude NH	3.5 1.142	2.6 1.6	10 ± 2
		4.0 1.117	0.6 2.3	6 ± 1
		4.5 1.116	0.5 1.6	10 ± 2
		5.0 1.124	1.3 0.53 38 ± 10
		5.5 1.132	2.0 0.11	440
		6.0 1.132	2.0 0.03	1791
		6.5 1.132	2.0 0.02	2977
	Drude tgNH 4.5 1.116	0.6 1.3	11 ± 3 8.6 ± 0.3

Table 6 .

 6 3 -Experimental (Mariani et al., 2017;[START_REF] Russina | Structure of a Binary Mixture of Ethylammonium Nitrate and Methanol[END_REF] and calculated equilibrium and dynamic properties of EAN at 298 K for different σ ON-N1 values.

		σ ON-N1	ρ	ρ dev D cat D ani	η eq
	Exp		1.210	3.86 7.80	32.7
	Drude NH	3.10	1.181 -2.36
		3.15	1.184 -2.15
		3.20	1.206 -0.28 0.05 0.08	2964
		3.25	1.213	0.32
		3.35	1.261	4.24
		3.50	1.267	4.77
		3.70	1.228	1.54
		3.75	1.217	0.61 4.15 6.92 52 ± 35
		3.80	1.205 -0.41 5.22 6.88 64 ± 19
		4.00	1.147 -5.21

Table 6 .

 6 4 -Equilibrium and dynamic properties of ChCl-EG predicted with different LJ parameters for the chloride anion at 298 K. ρ ρ dev D cat D ani D HBD η noneq

	Exp a	1.120	2.62	4.77	37
	CL&P	1.160 3.5	0.12 0.02	0.02	500 ± 105
	CL&P ref b 1.129 0.8	2.18 6.15	13.0	38 ± 2
	OPLS-AA c 1.130 0.9	1.73 3.62	5.61	45 ± 3

Table 6

 6 Units are: ρ/g cm -3 , ρ dev /%, (a, b, c)/Å.

	.6 -Experimental (Henderson et al., 2012) and calculated properties of ethylammonium nitrate
	crystal at 260 K.
		ρ	ρ err	a	b	c	α	β	γ
	Exp	1.367		39.6 46.1 39.9 90.0 112.7 90.0
	Drude NH	1.255 -8.2 39.5 50.0 39.4 90.0 110.3 90.1
	Drude tgNH 1.253 -8.3 39.6 50.0 39.4 90.0 110.6 90.0

  O (Ch + )...Cl - H OG (EG)...Cl - H 1 (Ch + )...Cl - H OG (EG)...O HG (EG) O (Ch + )...Cl - H U (U)...Cl - H 1 (Ch + )...Cl - H U (U)...O U (U) Figure 6.9 -Radial distribution functions of Cl H O , Cl H OG/U ,Cl H 1 , O HG/U H OG/U of ChCl EG at 298 K (left) and ChCl U at 353 K (right). The distributions of O H H O , O HG/U H O and O H H OG/U were not included in the plot due to absence of a significant signal.

		18				15
	g(r)	0 6 12	0	2	4	6 H 0 8 5 10 g(r)	0	2	4	H	6	8
					r / Å				r / Å

Table 6 .

 6 7 -Experimental[START_REF] Matsumoto | Thermal, Physical, and Electrochemical Properties of Li[END_REF] and calculated equilibrium and dynamic properties of[Li][FSI] [C 2 C 1 im][FSI] at 338 K.

		ρ	ρ dev D Li + D im + D ani	η eq	η non-eq
			x [Li][FSI] = 0.2
	Exp	1.460		11.1
	Drude NH	1.461 0.03 4.48 8.36 6.24	21 ± 3
	Drude tgNH 1.480 1.34 2.73 4.15 3.47	40 ± 3	33.3 ± 1.3
			x [Li][FSI] = 0.4
	Exp	1.536		21.8
	Drude NH	1.546 0.66 1.69 3.11 1.83	118 ± 61
	Drude tgNH 1.563 1.76 0.74 1.18 0.77 676 ± 417 143 ± 10

Table 7 .

 7 1 -System compositions given as numbers of molecules or ions in the simulation boxes.

	Component DES LiCl DES LiCl W
	Ch +	480	462
	Cl -	530	512
	U	960	924
	Li +	50	50
	W		275

  2.1-2.2 Å and a Cl. . . H O angle of 155-180°, depicted inAppendix (Fig. D.1). A weaker and less probable interaction is observed with the hydrogen atoms of urea (Fig.7.3b, grey curve) at 2.35 Å. The methyl and methylene hydrogens (green curve) are found at 2.9 Å from chloride and show broad and flat angular distributions in Fig.D.1 arising from a variety of available choline chlo-ride configurations, recently studied via first principles calculations[START_REF] Ashworth | Doubly ionic hydrogen bond interactions within the choline chloride-urea deep eutectic solvent[END_REF]. The first coordination shell includes 0.67 H O from cholinium (r min = 3.28 Å), 0.85 H w from water (r min = 3.02 Å) and 2.95 from H U urea (r min = 3.15 Å); the peak at 3.51 Å on the RDF of the H w of water was not included in the integration because it represents the second hydrogen atom from a water molecule, which is not interacting directly with chloride anion. Figure 7.3 -Radial distribution functions of selected hydrogen atoms around the chloride anion in the ChCl U LiCl W system. The H O and H 1 atoms correspond to the H atoms of hydroxyl and methyl groups of cholinium, and H U and H w are the H atoms of urea and water, respectively.

		12							4
			a				H O (Ch + ) H w (W)			b				H 1 (Ch + ) H U (U)
		9							3
	g(r)	6						g(r)	2
		3							1
		0	0	2	4	6	8	10	0	0	2	4	6	8	10
					r Cl -...X / Å							r Cl -...X / Å

Table 7 .

 7 2 -Experimental and calculated properties of DES-based electrolytes at 353 K.

	System	ρ	ρ dev D(Ch + ) D(Cl -) D(U) D(Li + ) D(W)	η
	DES LiCl W
	Exp	1.1662		8.1		15.2	6.1		22.5
	Sim	1.1402 -2.2	10.4	12.7	18.1	7.3	34.9	16 ± 2
	DES LiCl
	Sim	1.1645		3.4	4.9	7.4	3.1		55 ± 28

Table 7 .

 7 3 -Experimental and calculated reorientational times at 353 K.System τ(N(CH3 ) 3 ) τ(CH 2 N) τ(CH 2 O) τ(OH Ch ) τ(NH 2 ) τ(OH W ) τ(Li)

	DES LiCl W
	Exp	3.6	3.9	3.7		6.8	26.7
	Sim	4.9	3.3	3.8	8.4	7.6	5.0
	DES LiCl
	Sim	9.9	6.7	8.4	21.2	18.9
	Units are: τ/10 -11 s.

Table A .

 A 1 -Results of the fit of the ACF tail of p xy pressure component of [C 2 C 1 im][DCA] 303 K in the range 0.8-36.9 ps with the exponential functions.

	Function	τ	a	α	β	R 2	σ
	f 0 (t) f 1 (t) f 2 (t) f 3 (t) f 4 (t)	0.192 ± 0.009 6.9 ± 0.5	0.0110 ± 0.0003 0.0279 ± 0.0002 0.04 ± 0.01	< 0 0.9344 8.4 × 10 -4 5.0 × 10 -3 0.353 ± 0.005 0.9914 3.0 × 10 -4 0.2196 ± 0.0005 0.106 ± 0.002 0.9864 3.8 × 10 -4 0.75 ± 0.12 0.29 ± 0.03 0.9921 2.9 × 10 -4

  2 -The fitted parameters of eq. 2.65 applied to the running integrals of the p xy , p xz ,p yz pressure components of [C 2 C 1 im][BF 4 ] at 323 K valuated from a single trajectory (with and without fit constraints) and averaged over five independent trajectories. Units are: τ/ps, η/mPa s.Table A.3 -Viscosity values of ethylene glycol (EG) and ethylammonium nitrate (EAN) at 298 K atdifferent acceleration strengths, obtained using the periodic perturbation method.Units are: A/10 -5 Å fs -2 , η/mPa s. Table B.1 -Induction and dispersion contributions, and fraction of dispersion in non-Coulomb attraction, k ij , obtained with different SAPT levels and by the present predictive scheme, for the dimers of fragments considered in this work. Units are: E/kJ mol -1 . Table B.2 -Energy averages from MD trajectories corresponding to Lennard-Jones and electrostatic contributions for ionic liquids and isolated ions, used to calculate cohesive energies. FixQ: fixedcharge CL&P force field; Drude: CL&P model with polarisation added; SDrude: Drude model with scaled LJ by k SAPT2+ ij . Note the 0 values for BF - 4 because intramolecular interactions are excluded in the force field (there are no atoms separated by more than 2 bonds). Units are: E/kJ mol -1 .

	Dimer		E sSAPT0 tot	E sSAPT0 ind	E sSAPT0 disp	k sSAPT0 ij	E SAPT2+ tot	E SAPT2+ ind	E SAPT2+ disp	k SAPT2+ ij	pred k ij
						cation-anion
	C 2 C 1 im + • • • dca -C 2 C 1 im + • • • Ntf -2 C 2 C 1 im + • • • BF -4 E IL N	-358.6 -348.1 -355.2 E IL LJ N	-35.2 -29.8 -34.9 E IL elst N	-52.9 -51.9 -28.1 E +	0.60 0.64 0.45 E + LJ	-368.0 -362.1 -368.4 E + elst	E -	-40.8 -38.5 -42.7	E -LJ	-63.8 -72.3 -46.3 E -elst	0.61 0.65 0.52	0.68 0.55 0.50
	C 2 C 1 im + • • • MsO -	-415.8	-59.8 [C 4 C 1 im][BF 4 ] 343 K -36.5 0.38	-422.2		-72.4		-56.1	0.44	0.52
	C 2 C 1 im + • • • OAc -FixQ -153.7	-465.5 -63.5	-151.8 -283.9 306.2 -43.0	0.22 -1.3	-465.2 150.5	-177.3 22.8	0.0	-64.1 0.0	0.27	0.44
	C 2 C 1 im + • • • OTf -Drude -159.5	-362.8 -60.8	-44.9 -308.3 240.3 -36.2	0.45 -2.2	-375.9 125.5	-56.5 23.1	0.0	-54.7 0.0	0.49	0.45
	C 1 C 1 pyrr + • • • dca -SDrude -143.7	-338.6 -36.3	-33.0 -319.1 238.3 -34.1	0.51 -1.8	-349.6 125.3	-38.2 23.1	0.0	-42.9 0.0	0.53	0.52
	C 1 C 1 pyrr + • • • Ntf -2 C 1 C 1 pyrr + • • • BF -4 C 2 Py + • • • BF -4 FixQ -244.8 C 2 Py + • • • Ntf -2 Drude -254.9 N + 1111 • • • BF -4 SDrude -240.5	-332.7 -361.0 -358.2 -73.6 -342.5 -69.9 -358.1 -48.4	-40.2 -40.2 [C 2 C 1 im][dca] 303 K -38.8 0.49 -24.3 0.38 -35.1 -28.6 0.45 -304.7 240.3 -1.1 -37.0 -48.5 0.57 -348.8 244.2 -0.98 -37.0 -20.0 0.35 -360.6 244.2 -0.98	-342.9 -375.3 -373.6 141.7 -353.2 139.3 -0.61 -48.1 -49.1 -43.5 7.0 -0.50 -44.9 -0.56 -372.1 -45.6 139.3 -0.61 -0.56	-55.7 -41.3 -45.9 -16.8 -66.1 -26.3 -34.6 -26.3	0.54 0.46 0.51 0.60 0.43	0.50 0.40 0.49 0.56 0.38
	N + 1111 • • • Ntf -2	-334.4	-34.7 [C 4 C 1 im][Ntf 2 ] 323 K -36.6 0.51	-349.9		-44.1		-53.1	0.55	0.49
	FixQ C 2 C 1 im + • • • C 4 H 10 122.9 Drude 103.5 C 2 C 1 im + • • • C 6 H 14 SDrude 137.0	-119.5 -19.8 -116.4 -23.0 -80.8	-26.1 295.3 -9.2 -27.6 cation-neutral -1.5 150.5 0.75 -24.7 307.2 -10.3 -0.45 -83.5 234.2 -2.1 125.5 288.4 0.35 -10.8 -28.6 0.73 -29.2 -12.5 -87.2 234.0 -1.7 125.4 288.4 0.35	193.5 -33.2 152.3 -34.9 152.3	0.76 0.74	0.78 0.75
	C 2 C 1 im + • • • Bz	-49.3	-22.2 [C 4 C 1 pyrr][Ntf 2 ] 343 K -33.6 0.60	-47.7		-23.5		-37.8	0.62	0.74
	C 2 C 1 im + • • • Tol FixQ 317.5	-54.5 -91.5	-28.4 20.9 478.1 -45.4	0.61 23.4	-53.9 195.3 314.5 -26.6 -0.31	-44.7 194.2	0.63	0.75
	C 2 C 1 im + • • • DMSO Drude 295.0	-105.1 -89.0	-36.1 -38.7 413.1 -19.9	0.36 19.1	-88.2 159.1 294.7 -39.1 0.38	-28.4 153.1	0.42	0.57
	C 1 C 1 pyrr + • • • C 4 H 10 SDrude 338.0	-13.2 -44.7	-8.4 -41.1 406.3 -15.0	0.64 19.8	-17.6 159.3 294.7 -10.6 0.38	-21.8 153.1	0.67	0.69
	C 1 C 1 pyrr + • • • DME C 2 Py + • • • C 4 H 10 N + 1111 • • • C 4 H 10 -190.0 FixQ Drude -201.8	-47.5 -17.7 -12.7 -81.9 -77.4	-14.2 -7.8 [C 4 C 1 im][dca] 323 K -17.8 0.56 -23.6 0.75 -10.7 -14.4 0.57 -288.7 295.3 -1.5 -340.7 234.2 -2.1	-47.0 -21.3 -17.0 150.5 125.5	-16.5 -8.8 -12.1 8.2 -0.51 0.37 -0.56	-24.2 -27.3 -18.6 -16.7 -26.3	0.59 0.76 0.60	0.64 0.77 0.64
	SDrude	-182.5	-50.4	-354.5 234.0	anion-neutral -1.7 125.4	0.37	-0.56	-26.3
	Component C 4 H 10 • • • dca -C 6 H 14 • • • dca --14.1 FixQ C 4 H 10 • • • Ntf -2 Drude -28.1 C 6 H 14 • • • Ntf -2 SDrude -3.3	A -18.0 -25.1 -57.1 -17.2 -55.1 -26.5 -25.1	α [C 4 C 1 pyrr][dca] 323 K τ 1 -9.6 -17.2 0.64 -16.7 -21.8 0.57 -241.1 471.1 23.2 -9.2 -17.5 0.65 -286.6 400.2 18.4 -13.3 -25.2 0.65 -295.3 410.3 19.4	-25.6 -32.9 195.1 -25.0 159.2 -30.3 159.7	τ 2 -12.2 -18.4 8.2 -0.51 -10.5 0.37 -0.56 -10.4 0.37 -0.56	R 2 -26.8 -28.9 -16.7 -26.1 -26.3 -34.4 -26.3	σ 0.69 0.61 0.77 0.77	0.72 0.71 0.69 η 0.72
	1 trajectory p xy p xz p xz cons p yz 5 trajectories p xy ave C 4 H 10 • • • BF -4 C 6 H 14 • • • BF -4 FixQ -154.2 0.64 ± 0.02 0.9897 ± 0.0005 14.9 ± 0.6 0.9 ± 0.1 0.83 ± 0.02 4 ± 1 0.67 ± 0.03 0.984 ± 0.002 13.8 ± 0.9 1000.00 ± 0.03 0.9846 0.64 83 ± 6 0.9817 0.79 > 10 7 0.9892 0.53 81388 15.7 19.7 0.52 ± 0.02 0.970 ± 0.003 15 ± 1 583 ± 99 0.9872 0.55 17.0 0.67 ± 0.03 0.913 ± 0.005 10.3 ± 0.7 168 ± 9 0.9956 0.37 16.2 0.65 ± 0.02 0.914 ± 0.004 11.7 ± 0.7 210 ± 10 0.9967 0.37 18.9 -20.0 -16.6 -11.2 0.40 -28.8 -19.1 -20.2 0.51 0.49 -26.1 -20.8 -13.9 0.40 -36.4 -23.7 -24.7 0.51 [C 6 C 1 im][dca] 323 K 0.50 -91.2 -282.9 339.7 -0.84 153.3 8.2 -0.51 -16.7 C 4 H 10 • • • MsO --28.3 -20.6 -18.3 0.47 -37.2 -23.1 -29.3 0.56 0.59 Drude -167.7 -86.6 -340.4 270.1 -2.0 115.7 0.37 -0.56 -26.3 DMSO • • • Ntf -2 -53.5 -12.4 -23.3 0.65 -58.0 -15.0 -34.0 0.69 0.60 SDrude -147.6 -56.9 -357.2 271.3 -1.4 115.6 0.37 -0.56 -26.3 DMSO • • • OAc --97.8 -39.5 -21.2 0.35 -102.1 -49.3 -35.9 0.42 0.48 DME • • • FSI --25.8 -7.5 -16.5 0.69 -34.1 -9.8 -26.7 0.73 [C 2 C 1 im][Ntf 2 ] 323 K 0.66 Tol • • • MsO --38.0 -26.5 -19.7 0.43 -44.8 -30.1 -30.6 0.50 FixQ 82.1 -109.8 -38.2 249.1 -0.94 141.8 307.2 -0.45 193.5 0.48 AN • • • Ntf -2 -35.7 -6.0 -11.1 0.65 -38.9 -7.0 -17.1 0.71 Drude 63.2 -106.8 -93.0 247.9 -0.92 139.4 288.4 0.35 152.3 0.55 Bz • • • OAc --50.0 -41.4 -16.5 0.29 -55.0 -46.6 -27.5 0.37 0.44 SDrude 93.7 -74.8 -96.3 247.9 -0.92 139.4 288.4 0.35 152.3 p xz ave p yz ave 1.3 ± 0.5 0.81 ± 0.06 2 ± 1 62 ± 6 0.9644 1.27 17.5 Bz • • • OTf --10.0 -4.5 -5.9 0.57 -14.5 -5.5 -9.9 0.64 0.43 [C 6 C 1 im][Ntf 2 ] 343 K
	FixQ C 4 H 10 • • • C 4 H 10 181.0 Drude 160.9 C 4 H 10 • • • Tol SDrude 196.5	-125.3 -5.5 -122.2 -15.4 -84.4	-19.9 348.8 neutral-neutral -0.78 153.0 -0.60 -9.3 0.94 -8.2 314.5 -0.89 -0.31 -80.0 276.2 -2.1 115.5 294.7 0.38 -2.9 -27.8 0.91 -16.6 -5.0 -83.8 279.0 -1.3 115.4 294.7 0.38	194.2 -14.1 153.1 -43.1 153.1	0.94 0.90	1.00 1.00
	DMSO • • • DMSO	-33.9	-11.5 [C 2 C 1 im][BF 4 ] 323 K -20.4 0.64	-34.4		-13.1		-30.6	0.70	0.68
	DMSO • • • Tol FixQ -212.3	-25.2 -55.7	-10.3 -300.5 249.1 -25.4	0.71 -0.94	-26.5 141.8	-11.2 22.4	0.0	-33.9 0.0	0.75	0.81
	DMSO • • • C 4 H 10 Drude -216.7	-11.2 -53.5	-6.0 -319.3 247.9 -16.0	0.73 -0.92	-15.8 139.4	21.8	-6.3	0.0	-23.2 0.0	0.79	0.81
	DME • • • DME SDrude -206.7	-14.0 -35.4	-3.0 -328.7 247.9 -18.5	0.86 -0.92	-19.6 139.4	21.8	-3.4	0.0	-26.6 0.0	0.89	1.00
	DME • • • C 4 H 10 AN • • • C 4 H 10 AN • • • AN FixQ -106.6 Bz • • • Bz Drude -140.1 SDrude -124.4	System EG -9.7 -6.0 -26.6 -71.0 -277.9 354.7 A η noneq 0.005 8.9 ± 1.4 0.01 9.3 ± 0.7 -1.6 -15.9 -1.5 -9.0 [C 6 C 1 im][BF 4 ] 353 K System EAN 0.91 -14.2 A 0.01 28.7 ± 2.4 η noneq 0.02 32.8 ± 1.5 -1.7 -21.6 0.86 -8.6 -1.7 -11.9 -5.8 -10.9 0.65 -25.3 -7.2 -16.2 -0.50 153.3 23.9 0.0 0.0 -5.7 -1.1 -21.2 0.95 -8.1 -1.1 -21.6 -71.7 -308.7 281.7 -1.9 115.4 23.3 0.0 0.0 -45.3 -322.0 281.2 -1.4 115.4 23.3 0.0 0.0	0.93 0.87 0.69 0.95	1.00 0.77 0.63 1.00
				0.02 8.6 ± 0.3	0.05 23.4 ± 0.4
				0.025 8.3 ± 0.2
				0.05 7.9 ± 0.1
				0.1	6.4 ± 0.1

2.4 dynamic properties as a criterion of reliability of the force field

joined me for the kitchen talks about literally everything.

I spent good