Structural, hydrogen bonding and dipolar properties of alkyl imidazolium-based ionic liquids: a classical and first-principles molecular dynamics study

Ionic liquids (ILs) feature a tailorable and wide range of structural, chemical and electronic properties that make this class of materials suitable to a broad variety of forefront applications in next-generation electronics. Yet, their intrinsic complexity call for special attention and experimental probes have still limitations in unraveling the interactions occurring both in the bulk IL and at the interface with the solid substrates used to build the devices. This works provides an atomistic insight into these fundamental interactions by molecular modeling to complement the information still not accessible to experiments. In particular, we shed some light on the nature of the chemical bonding, structure, charge distribution and dipolar properties of a series of alkyl-imidazolium-based ILs by a synergy of classical and first-principles molecular dynamics simulations. Special emphasis is given to the crucial issue of the hydrogen bond network formation ability depending either on the nature of the anion or on the length of the alkyl chain of the cation. The hydrogen bond strength is a fundamental indicator of the cohesive and ordering features of the ILs and, in this respect, might be exploited to foster a different behaviour of the IL used a bulk medium or when used in electronic devices.


I. INTRODUCTION
Ionic liquids (ILs) bear great promise in green chemistry, environmental science and innovative electronic devices due to their unique properties, such as a tailorable structure, chemical variety, and their environmentally friendly features [2][3][4]57 . ILs are particularly promising for the realization of devices able to trigger interesting phenomena such as superconductivity [4][5][6] and quantum interference 7 in 2D semiconducting (nano)materials [8][9][10][11] . ILs consist of cations and anions and are liquids near room temperature, they have attracted broad attention from the scientific and industrial communities because of their sizeable structural chemistry and configurable properties. Due to their rich interactions and microstructures, ILs may work as solvents, materials, and catalysts, rendering them broadly useful in green chemistry and novel electronic devices. Alkylimidazolium ionic liquids, especially ILs based on 1-alkyl-3-methylimidazolium cation, have been among the most studied ILs because of their remarkable properties and the simplicity of synthetic methods for producing them. These ILs features strong coupling interactions between the electrostatic forces and hydrogen bonds with complex nano-and micro-structural ordering [12][13][14] . In particular, the peculiar properties and functions of ILs are dominated by a unique class of hydrogen bonds. [15][16][17] Such complexity facilitate the formation of complex nanostructures (including ionic pairs, aggregates, or even ionic clusters) thus, the corresponding efa) Electronic mail: mauro.boero@ipcms.unistra.fr b) Electronic mail: guido.ori@ipcms.unistra.fr fects on the performance of ILs in typical applications are significant. 14 The complex interplay between hydrogen bonds (H-bond) and dipole properties of ILs is straightly dependent on the features of the IL used. Beside their tunability in terms of chemical compositions and physicochemical properties, these features have been identified as affecting ILs performance whether used as bulk solvent phases for catalysis or as in gating media in optoelectronic transistors [18][19][20][21][22][23] . Atomistic modeling, especially in terms of classical molecular dynamics (MD) simulations, has been very useful in the last decade in improving our understating of ILs behaviour by complementing experimental data [24][25][26][27] . Several examples of MD studies oriented to the comprehension of H-bonds features in ILs are available nowadays. [28][29][30][31][32] Ab initio MD has also been employed to obtain further insight into the structure, Hbonds and electronic properties of a few ILs [33][34][35][36][37][38][39][40] . However these works have been focused on individual or short series of ILs. The purpose of this work is to explore via molecular dynamics (MD) methods structural, hydrogen bonding and dipole properties of a large series of alkylimidazolium-based IL commonly used in optoelectronic devices. The specific family of ILs on which we focus is constituted by imidazolium-based cations with different lengths of the alkyl chain and the mostly experimentally exploited inorganic anions, as sketched in Figure 1.
Our computational study makes use of classical molecular dynamics (MD) to simulate, with the same consistent approach, a large numbers of IL systems and to extract their structural properties and related H-bond network features. For a selected number of IL systems, first-principles molecular dynamics (FPMD) simulations are used to enrich the information by the inclusion of the electronic structure to refine our understanding of their interaction without relying on ad-hoc parameterized force fields. This allow for a more precise understanding of the chemical bonding and the electrostatics regulating the H-bond network formation and dipolar interactions in comparison with the results obtained by classical MD. The ability of our FPMD approach in simulating these ILs at finite temperature has already been assessed 41 , showing that this is a reliable computational protocol for the purpose of our present investigation. The novelty of the study presented here stems from the relatively large number of systems considered (four cations of different alkyl chain length combined with four different anions), the accurate comparison of their peculiar properties and the complementary electronic properties quantified in terms of dipole moments, accessible by FPMD, responsible for the type and topology of H-bond network realized by these compounds. The relative strength of the H-bond depends on the chemical nature and molecular structure of cations and anions and is a fundamental indicator of the cohesive features of the ILs.

A. Classical molecular dynamics
The ILs on which this work is focused are made of imidazolium-based cations, specifically 1-alkyl-3methylimidazolium cations, combined with different inorganic anions having a general formula C x MIM + -A − , where x = 1, 2, 4 and 6 (also termed DMIM, EMIM, BMIM, and HMIM, respectively) and the anion A − is chlorine Cl, or tetrafluoro-borate BF 4 , or hexafluoro-phosphate PF 6 or bis(trifluoromethane)-sulfonimide TFSI. On a general ground, a nonpolarizable force fields (FF) is sufficient to reproduce with appreciable accuracy structure and local bonding of the targeted ILs 29,42-44 . In the present work, we use the FF originally proposed and developed by Canongia Lopez and coworkers 24 with the parametrization of Koddermann et al. 45 to better account for the dynamical properties of the bulk ILs simulated. This specific FF is fitted on ab initio calculations and allows to simulate a wide variety of ILs chemical structures including the specific ones targeted in this work.
Periodic boundary conditions are applied to the bulk IL simulation cell and the dispersive interactions are neglected beyond a cutoff of 15Å. A standard Ewald summation is used in the calculation of electrostatic interaction to avoid finite size effects, with the Ewald parameters selected to provide an accuracy of 10 −5 atomic units (au). For each IL, we constructed two simulation cells of 20 and 200 IL pairs. This size difference will be clarified in the next paragraph and in the ongoing discussion. The larger systems accounting for 200 IL pairs are sufficiently large to ensure a statistically significant average of the structural properties (See Table S1 in the Supplementary Material). The initial configuration has been realized as a random accommodation of IL pairs inside a cubic box by means of the packmol package 46 and the initial volume of this box was selected based on the experimental density of the IL. The bulk model of the IL was then produced following a standard melt-quench procedure heating the systems up to ∼ 1000 K for 2.5 ns then by cooling these same systems to room temperature (∼ 300 K) in a sequence of steps in which a temperature decrease of 100 K was imposed and on which the system was allowed to re-equilibrate for 1 ns in the canonical ensemble (NVT). ns in the canonical ensemble (NVT). The initial 1000 K NVT stage lasted for 2.5 ns to ensure the lose of memory of the intial configuration constructed as reported above. After the NVT thermal cycle, a simulation stage in the isobaric-isothermal ensemble (NPT) at T = 300 K and P = 1 atm was performed for 1.5 ns to get well equilibrated models and to check the density of the simulated systems with respect to the experimental values. The outcome of these simulations confirmed that the theoretical density obtained for our model ILs has small deviations from the experimental value (with the exception of DMIM-Cl (∼5%), the majority showed deviaitons <3%). These preparatory and relatively longlasting simulations provided the actual systems simulated for about 10 ns in the NVT ensemble. This is the production run stage in which statistics was generated and used to fully characterize the various ILs. More precisely, the data analyzed and presented in the following paragraphs refer to the final 5 ns in the NVT ensemble where trajectories were sampled each 1000 steps. Temperature and pressure were controlled using a Nosé-Hoover thermostat [47][48][49] . The simulations were performed with the DLPOLY package 50 , and a standard Verlet numerical integration scheme was used with a time step of 1.0 fs.

B. First-principles molecular dynamics
To account for the effects of an explicit inclusion of the electronic structure, we performed FPMD simulations of a few IL model systems. To this aim, to make the system tractable within this type of approach and to mitigate the computational workload, we used the smaller IL systems mentioned in the former paragraph. Our FPMD is based on the density functional theory (DFT) of Kohn and Sham 51 (DFT) and the dynamical simulations were performed according to the Car-Parrinello method 52 as implemented in the developer's version of the CPMD 53 code. The valence-core interaction was described by norm-conserving numerical pseudopotential of the Troullier-Martins 54,55 type, and the exchange and correlation interactions are expressed according to the functional proposed by Becke 56 and Lee, Yang, and Parr 58 (BLYP), respectively. Long-range van der Waals dispersion interactions, not included in the BLYP functional, are added according to the maximally localized Wannier functions 59 and centers (WFCs) approach 60 boosted by a propagation scheme 42 to reduce the computational burden. Valence electrons are treated explicitly and their wavefunctions are expanded on a plane-wave (PW) basis set with a cut-off of 100 Ry (1360.57 eV). This value was tuned based on the convergence of the stress tensor (pressure) as a function of the number of PWs in the preparatory stage before moving to actual production run (See Figure S1 in the Supplementary Material). We remark that this relatively large cut-off is crucial for an accurate description of the electron density around several chemical species (O, N, S, F, Cl) composing the ILs. Canonical NVT simulations were done by controlling the temperature via a Nosé-Hover 47-49 thermostat chain. 61 For the numerical integration of the Car-Parrinello equations of motion a fictitious electron mass of 400.0 au and a time step of 3.0 au (0.0726 fs) ensured good numerical control of the constants of motion. The reliability of our CPMD computational protocol for the simulation of this class of systems has already been reported elsewhere. 41 The CPMD simulated ILs ([DMIM][Cl, BF 4 , PF 6 and TFSI] and [EMIM][BF 4 ]) started from liquid phase models amounting to 20 pairs, which translates into a number of atoms between 340 and 620, depending on the specific composition. A pre-equilibration stage within our classical MD scheme and following the thermal cycle explained above was the procedure adopted to generate the initial configurations of the present set of FPMD simulations. The stress tensor analysis done during NVT dynamical simulations at different cell sizes allowed to assess the actual simulation cell to be used for the production and the resulting theoretical density of the ILs turned out in fair agreement (< 4%) with the experimental value. The system displays just a minimal residual pressure (< ± 4 GPa), (see Figure S1 of the Supplementary Material), compatible with the DFT accuracy. 62

A. Structural properties
A first issue on which we focused was the structure of the alkyl-imidazolium ILs quantified in terms of cationanion pair correlation functions (PCFs) g ij (r) using the standard definition where N is the number of molecules simulated inside the cell of volume V and r ij is the distance between the center of mass (COM) of each cation and anion constituting the IL. The result of these COM PCFs for the different ILs considered here are reported in Figure 2. We remind that no size-effect have been found comparing 20 and 200 ILs pairs systems with classical MD when analysing the first coordination shells around the ions (see Figure S2 in the Supplementary Material). These results give an insightful global picture of the structures as a function of both the cation alkyl chain and the type of inorganic anion. For the ILs in which Cl − is the anion, in the first coordination shell of the cation-anion pairs, the shortest alkyl chain (DMIM) shows a broad peak centered at 4.5Å, followed by a smaller peak at ∼ 5.8Å. Then, by increasing the length of the alkyl chain from one to six carbon atoms, the first peak displays an increasing broadening until the two peaks nearly merge in a unique feature for the longer HMIM, as an effect of the larger steric volume characterizing this cation. For the BF 4 series, a similar trend is observed but overall the intensity of the peaks is reduced with respect to the PCFs of the former ILs carrying as an anion Cl − . The PCFs of DMIM-and EMIM-BF 4 show a double peak with the two maxima centered at 4.5 and 5.2Å, a clear indication of two preferred local coordination of the anion with respect of the orientation of the alkyl-imidazolium cation. In the case of BMIM, the double peak is replaced by a broad first peak is centered at 4.4Å together with minor second peak at 6.2Å. Finally, the PF 6 and TFSI compounds show analogous trends with a further reduction in the intensity with respect to both Cl − and BF 4 anions. Namely, a double peak present in the case of the short alkyl chains becomes a unique broad peak upon increasing of the chain length, as seen in the former cases. FPMD simulations performed on the same ILs models (DMIM) show PCFs trends similar to those obtained by the present classical MD approach in terms of position of the first peak. A major difference is that instead of a double peak, DFT-based dynamical simulations display a single first peak but more structured than the one obtained by classical MD. An additional insightful parameter characterizing the structure of the ILs is the coordination number n c . This can be obtained by integrating the PCF, g ij (r), as where ρ 0 and r c are the atomic density of each IL pairs and the position of the first minimum minimum of the PCF, respectively. The values of n c obtained from both classical MD and FPMD are reported in the table (I).
Using the cutoffs reported in Table ??, the calculated total average coordination numbers for the cation-anion in the first coordination shell fall in the range from 6 to 10. The size and steric hindrance of the anion, which increases as Cl < BF 4 < PF 6 < TFSI, mostly determine the average value of the coordination number, which turns out to be n c = 6.69, 7.46, 7.  This provides and additional validation of the FF used here.

B. Hydrogen bonding
We inspected the details of the H-bond network characterizing the various ILs targeted in this work to get a clearer picture of their nature and influence on the structure of these systems. A fundamental driving force determining the molecular structuring of IL is the electrostatic interaction occurring between IL pairs, responsible for the realization of the H-bond network of the system and, thus, of the cohesive properties of the cation and anion moieties. This general observation is supported by the fact that along the evolution observed in the sampled trajectories, the distributions of the anions (Cl − , BF − 4 , PF − 6 and TFSI) around cations is mainly governed by electrostatics and H-bond features, as shown by Figure  3, where the PCF of the different H atoms constituting the EMIM cation and the F atoms of the BF − 4 anion are considered. Indeed, these are the two atoms that form the inter-molecular H-bonds in this specific compound and serve as a practical example of this type of interactions in an IL.
By examining the PCF of Figure 3, we can see that the sharp peak realized by H 2 , which is the H atom shared by the two N atoms of the imidazolium ring, in the g H2F (r) distribution is located at 2.05Å. This is a clear indication of a strong preference of F atoms for binding to this site. Such a short distance is also a fingerprint of a tight Hbond. ? ? The other H atoms of the cation are characterized by corresponding PCFs peaked at larger distances, namely 2.5-2.7Å and reduced intensities with respect to g H2F (r). These are still active H-bonds, concurring to the overall structure of the IL ? , but being somehow weaker, they can be broken and reformed in a dynamical way as observed in the case of [EMIM] [TFSI]. 41 Analogous results for the other ILs targeted in this study are reported in Figure S3    this case, corresponding to the shortest alkyl chain considered here, the differences in these distributions can be clearly ascribed to the nature of the anion. We remark generally well localized CDF peaks, with the strongest localization corresponding to the BF − 4 for a distance of about 2Åand an angle of 150 • . Broader distribution characterize, instead, the cases of Cl − and TFSI in terms of angles which span the intervals 100 • -160 • and 85 • -160 • respectively. This is not entirely surprising. Indeed, the small Cl − anion is more mobile and can find stable positions in a slightly wider angular range than BF − 4 . On the other hand, both the larger PF − 6 and the flexible TFSI molecule, the latter having also different conformers 41 , require more space to be accomodated beside DMIM cations and this accounts for a more dispersed CDF distribution.

C. Dipole moments
In the case of FPMD simulations, the availability of the electronic structure allows for a straightforward calculation of the dipole moment. To this aim, we make use of the maximally localized Wannier functions and centers (WFCs) 59 , consisting in a unitary transformation of the Kohn-Sham orbitals. The WFCs provide a useful tool to analyze the chemical bonds and a shorthand visualization of the electronic structure. An example of WFCs representation for the different cations and anions of the ILs targeted here is reported in Figure S4 of the Supplementary Material. The WFCs allow also for an easy calculation of the local dipole moment µ of each cation-anion pair constituting the IL in terms of a practical point-like charges sum, namely where N is the number of atoms in a IL ion, q i is the valence (in the case of pseudopotentials approaches) atomic charge of atom i, R i are the coordinates of atom i, n W F C is number of WFCs and f s is the occupation of the s th Wannier orbital, namely f s = 2 in a spinrestricted calculation and f s = 1 in the spin-unrestricted case, and r s are the cartesian coordinates of each WFC belonging to an IL cation or anion. The analysis of the value of the dipole moment µ = | µ| was performed on the basis of MLWFs computed over YYY uncorrelated configurations for each IL system sampled along the dynamical trajectories obtained at 300 K for . We recall that the only gauge invariant dipole having a physical meaning in charge neutrality condition of the simulation cell is the one between the cation and the anion. ? In our analysis, we computed the dipole moments between each cation-anion pair present in our simulation cells and the statistics is extracted based on both the number of pairs and the different configurations sampled along the trajectories. The result is shown in Figure 7.
The , where the average dipole moment takes the value of 33.0 D, consistently with the more bulky anion structure. A general picture that can be extracted is the relative stability and insensitivity of the local dipole moment to the type of anion (Cl, BF 4 , PF 6 ), provided that such an anion has not a complex molecular structure, for the same IL cation (DMIM). An analogous minimal effect has been found for different alkyl chains in the cases of DMIM  Figure 7 is mainly due to the more complex molecular structure of the [TFSI] with respect to the other smaller and roughly spherically isotropic anions (Cl, BF 4 , PF 6 ). The bigger TFSI molecule has a larger steric encumbrance and is accommodated at larger distances from the cation, as shown in the analysis of the PCFs of Figure 1. As a consequence, the relative distance between the center of mass of the positive and negative charges increases, thus accounting for an increased dipole moment, in line with the results obtained for [EMIM]-[TFSI]. 41 As expected, the presence of an extended H-bond network induces a significant enhancement of the dipole moment of a corresponding dimer in the gas phase (14. In Figure 7, we report in the inset also the values of µ as a function of the interionic distance of the cation-anion pair extending up to the second coordination shell (0-15 A). A double-peak feature characterizes the IL in which the anions are Cl − , BF − 4 and PF − 6 , with the two maxima centered at ∼ 5 and ∼ 9Å. Since the dipole moment of an ionic liquid is experimentally difficult to measure, this data can provide a useful quantitative analysis to complement experimental probes.

III. CONCLUSIONS
We have inspected a series of alkyl-imidazolium-based ILs, targeted in forefront electronic applications, with the scope of unraveling their atomistic structure, specific interactions and microscopic features. By resorting to a combination of classical and first-principles dynamical simulations, we have shown how the composition of the ILs influences the structure, the nature and topology of the H-bond network and, in turn, the electronic properties, quantified in terms of dipole moment distribution of the different ILs. Our study makes use of an identical simulation protocol, based on both classical model potentials and DFT-based first-principles molecular dynamics, of a series of ILs on which forefront research is focused. This comparative study, done within an identical and consistent paradigm is prone to give a useful guidelines in the selection of the best suited IL systems according to the specific application targeted where the ability of H-bond may play a significant role.

SUPPLEMENTARY MATERIAL
The Supplementary Material include data and supporting figures quoted in the text.

DATA AVAILABILITY
Representative configurations of the different IL obtained by classical MD and FPMD raw trajectory data are available at the European Center of Excellence Novel Materials Discovery (CoE-NOMAD) repository. 63