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

Iréné Amiehe Essomba  Université de Strasbourg, CNRS, Institut de Physique et
Chimie des Matériaux de Strasbourg, UMR 7504, F-67034 Strasbourg, France
Mauro Boero<EMAIL_ADDRESS>Université de Strasbourg, CNRS,
Institut de Physique et Chimie des Matériaux de Strasbourg, UMR 7504, F-67034
Strasbourg, France  Kerstin Falk  Fraunhofer IWM, Freiburg, Germany  Guido Ori
<EMAIL_ADDRESS>Université de Strasbourg, CNRS, Institut de
Physique et Chimie des Matériaux de Strasbourg, UMR 7504, F-67034 Strasbourg,
France

###### Abstract

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 featuresWang2021 ; Vioux2017 ; Welton2004 ; ye2012 . ILs are
particularly promising for the realization of devices able to trigger
interesting phenomena such as superconductivityMorpurgo2015 ; Morpurgo2016 ;
ye2012  and quantum interferenceJia2018  in 2D semiconducting
(nano)materialsOrgiu2016 ; PEI2022100159 ; adfm.201909736 ; Santato2020 . 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 orderingWelton2014 ; Canongia2006 ; Wang2022 . In particular, the
peculiar properties and functions of ILs are dominated by a unique class of
hydrogen bonds.Dong2021 ; Dereka2021 ; Hunt2015  Such complexity facilitate
the formation of complex nanostructures (including ionic pairs, aggregates, or
even ionic clusters) thus, the corresponding effects on the performance of ILs
in typical applications are significant.Wang2022  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 transistorsWilkes2004 ;
Nguyen2022 ; Matsumoto2017 ; Yadav2021 ; Park2021 ; Velpula . 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 dataCanongia2004 ;
Coasne2011 ; Ori2014 ; Ori2015 . Several examples of MD studies oriented to
the comprehension of H-bonds features in ILs are available
nowadays.Morales2011 ; Agrawal2018 ; Niazi2013 ; Bhattacharjee2020 ;
Andreussi2012  Ab initio MD has also been employed to obtain further insight
into the structure, H-bonds and electronic properties of a few ILsBuhl2005 ;
DelPopolo2005 ; Grimme2012 ; Hunt2006 ; Cremer2010 ; Kirchner2015 ;
Kirchner2010 ; Fayer2017 . 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 alkyl-imidazolium-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.

Figure 1: (Color online) A) Chemical structure of the ionic liquid cations and
anions simulated in this work: as cations 1-methyl(DMIM), 1-ethyl (EMIM),
1-butyl (BMIM) and 1-hexyl-3-methylimidazolium (HMIM) and, as anions chlorine
(Cl), tetrafluoro-borate (BF4), hexafluoro-phosphate (PF6) and
bis(trifluoromethylsulfonyl)imide (TFSI). B) Example of atom indexing for
1-ethyl-3-methylimidazolium [EMIM] cation. C) Typical molecular configurations
of 200 IL pairs simulated by MD. The blue spheres are the EMIM whereas green,
orange , yellow and purple correspond to Cl, BF4, PF6 and TFSI anions.

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 assessedKana22 ,
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.

###  I.1 Classical molecular dynamics

The ILs on which this work is focused are made of imidazolium–based cations,
specifically 1-alkyl-3-methylimidazolium cations, combined with different
inorganic anions having a general formula CxMIM+ -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 BF4, or hexafluoro-phosphate PF6 or
bis(trifluoromethane)-sulfonimide TFSI. On a general ground, a non-polarizable
force fields (FF) is sufficient to reproduce with appreciable accuracy
structure and local bonding of the targeted ILs Andreussi2012 ; Boero2015 ;
Kirchner2018 ; Holm2018 . In the present work, we use the FF originally
proposed and developed by Canongia Lopez and coworkers Canongia2004  with the
parametrization of Koddermann et al. Koddermann2013  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 packmol  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 $\sim
1000$ K for 2.5 ns then by cooling these same systems to room temperature
($\sim 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 ${\rm T}=300$ K
and ${\rm 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 ($\sim$5%), the majority
showed deviaitons $<$3%). These preparatory and relatively long-lasting
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 Nose1 ; Nose2 ;
Hoover . The simulations were performed with the DLPOLY package dl_poly , and
a standard Verlet numerical integration scheme was used with a time step of
1.0 fs.

###  I.2 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 KS1965  (DFT) and the dynamical simulations were performed
according to the Car–Parrinello methodCP1985  as implemented in the
developer’s version of the CPMDCPMDcode  code. The valence-core interaction
was described by norm-conserving numerical pseudopotential of the Troullier-
MartinsTMpseudo1 ; TMpseudo2  type, and the exchange and correlation
interactions are expressed according to the functional proposed by
Beckebecke1988  and Lee, Yang, and Parrlee-yang-parr1988  (BLYP),
respectively. Long-range van der Waals dispersion interactions, not included
in the BLYP functional, are added according to the maximally localized Wannier
functionsMarzari2012  and centers (WFCs) approachSilvestrelli2008  boosted by
a propagation schemeBoero2015  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é–HoverNose1 ; Nose2 ; Hoover  thermostat chain.nose1992  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.Kana22  The CPMD simulated ILs ([DMIM][Cl, BF4, PF6
and TFSI] and [EMIM][BF4]) 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 ($<\pm$ 4 GPa), (see Figure S1 of
the Supplementary Material), compatible with the DFT accuracy.DFTacc
Trajectories were accumulated in the NVT production run at 300 K for 39 ps for
[DMIM][Cl], 34 ps for [DMIM][BF4], 27 ps for [DMIM][PF6]), 15 ps for
[DMIM][TFSI] and 20 ps for [EMIM][BF4]. 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) repositoryNOMAD .

##  II Results and discussion

###  II.1 Structural properties

A first issue on which we focused was the structure of the alkyl-imidazolium
ILs quantified in terms of cation-anion pair correlation functions (PCFs)
${\rm g}_{ij}(r)$ using the standard definition

$g(r)=\frac{V}{4\pi Nr^{2}}\sum_{i}\sum_{i\neq
j}\left<\delta({r}-{r_{ij}})\right>$ (1)

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.

Figure 2: (Color online) Cation-anion center of mass pair correlation
functions (gC/A(r)) for the ILs simulated by classical MD at 300K (200 ILs
pairs). In each plot, the insets show the gC/A(r) obtained by FPMD simulations
(20 IL pairs).

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 $\sim 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 BF4 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-
BF4 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 PF6 and TFSI compounds show
analogous trends with a further reduction in the intensity with respect to
both Cl- and BF4 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

$\centering n_{c}=4\pi\rho_{0}\int_{0}^{r_{c}}r^{2}g(r)dr\@add@centering$ (2)

where $\rho_{0}$ and rc 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
(1).

Table 1:  Cation-anion coordination numbers $n_{c}$ obtained from classical MD and FPMD simulations at 300K.222 Minimum of the $g_{ij}(r)$ after the first peak used as cutoff for the calculation of $n_{c}$ (in parenthesis the values used for FPMD): C1,2,4,6MIM-Cl, 6.46(6.43) Å, 6.86 Å, 7.26 Å, 7.86 Å; C1,2,4,6MIM-BF4, 7.25(7.32) Å, 7.64(7.52) Å, 8.05 Å, 8.06 Å; C1,2,4,6MIM-PF6, 7.68(7.46) Å, 8.15 Å, 8.48 Å, 7.26 Å; C1,2,4,6MIM-TFSI, 8.98(8.98) Å, 9.17 Å, 9.58 Å, 9.98 Å. |  | Cx-Cl | Cx-BF4 | Cx-PF6 | Cx-TFSI
---|---|---|---|---|---
DMIM | MD | 6.69(3) | 7.46(3) | 7.92(3) | 10.11(2)
| FPMD | 7.01(4) | 7.44(2) | 7.12(3) | 10.02(2)
EMIM | MD | 6.99(2) | 7.48(3) | 7.99(2) | 7.51(3)
| FPMD | – | 7.32(3) | – | –
BMIM | MD | 6.10(3) | 6.91(2) | 7.17(4) | 7.30(2)
HMIM | MD | 6.77(2) | 5.93(3) | 6.57(3) | 6.36(2)

Using the cutoffs reported in Table 1, 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 $<$ BF4 $<$ PF6 $<$ TFSI, mostly determine the average value
of the coordination number, which turns out to be $n_{c}=$ 6.69, 7.46, 7.92
and 10.11 for DMIM-Cl,-BF4, PF6 and TFSI, in this order. The cation-anion
coordination number for Cl-based IL does not show any particular trend whereas
BF${}_{4}^{-}$, PF${}_{6}^{-}$ and TFSI show a decrease of $n_{c}$ upon
increasing the length of the alkyl chain. A figure of merit worthy of note is
the fact that the values of $n_{c}$ estimated from the classical MD show a
remarkable agreement with those obtained by FPMD. This provides and additional
validation of the FF used here.

###  II.2 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.

Figure 3: (Color online) Pair correlation functions gij(r) of the different H
atoms of [EMIM] with the F atoms of [BF4] obtained by classical MD at 300 K.
The atom numbering is the standard one given in the literature and reported
explicitly in Figure S2 of the Supplementary Material. Inset: a zoom-in in the
1-6 Å region.

By examining the PCF of Figure 3, we can see that the sharp peak realized by
H2, which is the H atom shared by the two N atoms of the imidazolium ring, in
the g${}_{H_{2}F}$(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 H-bond.HbondILs; HbondILs2 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${}_{H_{2}F}$(r). These are still active H-bonds, concurring to the overall
structure of the IL HbondILs2, but being somehow weaker, they can be broken
and reformed in a dynamical way as observed in the case of [EMIM][TFSI].Kana22
Analogous results for the other ILs targeted in this study are reported in
Figure S3 of the Supplementary Material. However, with respect to [EMIM][BF4]
some differences can be pointed out. In particular, [EMIM][Cl] PCFs shows that
together with H2 atom also the other H atoms within the imidazolium ring as
well as the H atoms of the first C atoms connected to the ring (C6 and C7)
show clear sharp peaks. These H-bond interactions are found at larger distance
with respect to to the case of F in [EMIM][BF${}_{4}7$], being at 2.65 Å
(H2,H4,5) and 2.85 Å (H6,H7) for respectively. In the case of [EMIM][PF6], the
interaction between F atoms and H2 atoms is found to be again the strongest
with respect to the other H atoms of the cation (2.25 Å versus 2.75–2.95 Å),
however the difference in intensity of the different H is much less pronounced
with respect to [EMIM][BF4]. In the case of [EMIM][TFSI], the strongest H-bond
interaction is found to be the one between H atoms of the cation and the O
atoms of the sulphonyl groups of TFSI. In this case, the H2-O PCF first peak
is found to be the most intense however all the Hx-O show a H-bond interaction
at a similar distance (2.66–2.7Å).
The effect of the alkyl chain length is summarized in Figure 4. A feature that
emerges is a structuring effect of the H-bond promoted by H2 with the increase
of the number of C atoms composing the alkyl chain, from DMIM to the longer
HMIM.

Figure 4: (Color online) Pair correlation functions gij(r) of the H atom
(labeled as H2) shared by the two N atoms of the IL cation and the Cl, F and O
atoms atoms of the anion. The data is based on the MD simulations at 300K and
is plotted for each anion as a function of the length of the cation alkyl
chain. Inset: data obtained by FPMD simulations.

This alkyl chain length effect is present in all the four ILs differing in the
anion moiety. The main difference is a decrease in the intensity of the first
peak according to the trend Cl$>$ BF${}_{4}>$ PF${}_{6}>$ TFSI. The comparison
with the FPMD results for the DMIM and EMIM systems, confirms the trend
evidenced by classical MD but reveals also the presence of slightly shorter
H2-Cl and H2-O distances in the cases of DMIM-Cl and DMIM-TFSI, respectively.
In Figure 5 we show the combined distribution function (CDF), in which the
horizontal axes represent the distance of the hydrogen atom H2 of the cation
form the F site of the BF${}_{4}^{-}$ anion participating to the formation of
H-bond. The vertical axes of each panel, instead, show the angle defined by
the vector oriented from H2 to the C2 atom of the cation and the vector
connecting H2 and F as sketched in Figure 5.

Figure 5: (Color online) Top: Combined distribution function showing the
hydrogen bond geometry between the atom H2 of the IL cation and the F atoms of
of the BF${}_{4}^{-}$ anion as a function of the length of the cation alkyl
chain (DMIM, EMIM, BMIM and HMIM) obtained by classical MD at 300 K. Bottom:
Left, Combined distribution function showing the hydrogen bond obtained by
FPMD; Right, Labeling of the main atoms quoted in the main text and vectors
defining the angle for which distribution density plots are shown in the
manuscript.

In the angle analysis presented in Figure 5 value of 180∘ indicates that the
atoms C2, H2 and F are aligned, thus the H-bond results perfectly linear. An
intense maximum arises in the region around 2 Å for a range of value between
125∘ and 180∘ for the IL pair on which we focus. This originates from the
H-bond formed by the H2 atom of the cation and and the acceptor F atom of the
BF4. By increasing the length of the alkyl chain, a narrowing of the
dispersion of H-bonding appears in the region around 2 Å and angles of
125∘-180∘. This is in perfect agreement with the features displayed by the
PCFs in Figures 3 and 4. The lower panels of Figure 5 shows the CDF data
obtained from FPMD data for DMIM-BF4 and EMIM-BF4 systems. The FPMD data show
a position of the strong H-bond at values similar to those computed from the
classical MD simulations, but the distributions are less localized and
affected by a larger broadening. This can be ascribed both to the oscillations
of the electronic structure and to the reduced statistics of FPMD simulations
in terms of system size and simulation time scales.

Figure 6: (Color online) Combined distribution function showing the hydrogen
bond geometry between the atom H2 of the DMIM and the Cl, F and O atoms of the
anions obtained by classical MD at 300 K.

Figure 6 reports analogous two-dimensional maps of the CDFs in the case of
different anions at fixed (DMIM) cation. Note that for BF${}_{4}^{-}$ and
PF${}_{6}^{-}$ the accpetor anion atom is F whereas for TFSI is the sulphonyl
O atom. In 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 conformersKana22 , require more space to be accomodated beside DMIM
cations and this accounts for a more dispersed CDF distribution.

###  II.3 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)
Marzari2012 , 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 ${\bf\mu}$ of each cation-
anion pair constituting the IL in terms of a practical point-like charges sum,
namely

$\centering{\vec{\mu}}=\sum^{N}_{i=1}q_{i}{\bf
R}_{i}-\sum^{n_{WFC}}_{s=1}f_{s}{\bf r}_{s}\@add@centering$ (3)

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$, ${\bf R}_{i}$
are the coordinates of atom $i$, $n_{WFC}$ is number of WFCs and $f_{s}$ is
the occupation of the $s^{\rm th}$ Wannier orbital, namely $f_{s}=2$ in a
spin-restricted calculation and $f_{s}=1$ in the spin-unrestricted case, and
${\bf 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
$\mu=|{\vec{\mu}}|$ 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 [DMIM]-[Cl], -[BF4], -[PF6], -[TFSI] and
[EMIM][BF${}_{4}^{-}$]. 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.Kirchner2011 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.

Figure 7: (Color online) Averaged dipole moment distributions of the IL pairs
inside the different ionic liquid. Inset: Averaged dipole moment distribution
plotted as a function of the distance of two counter ions.

The distribution of the dipole moment show an average value of 27.0 D with a
full width at half the maximum (FWHM) of about 10 D for both
[DMIM]-[BF${}_{4}^{-}$] and [DMIM]-[PF${}_{6}^{-}$]. In the case of
[DMIM]-[Cl-], while the FWHM does not undergo any change, the distribution
shifts toward slightly higher values with a peak centerd at 28.0 D. The most
remarkable displacement of the distribution was found for the case of
[DMIM]-[TFSI], 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, BF4, PF6), 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
and EMIM when BF4 is the cation. The larger dipole moment found in the case of
[DMIM]-[TFSI], as opposed to the other cases presented in 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, BF4, PF6). 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].Kana22

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.3 D, 13.6 D, 13.9 D, 14.1 D, 17.1 D, for [DMIM]–[Cl], [DMIM]–[BF4],
[DMIM]–[PF6], [DMIM]–[TFSI] and [EMIM]–[BF4], respectively). This enhancement,
typical of any polar liquid, is again an indicator of the presence of H-bonds
and a fingerprint of the crucial electrostatic interactions between cation and
anionKana22 .

In Figure 7, we report in the inset also the values of $\mu$ as a function of
the interionic distance of the cation-anion pair extending up to the second
coordination shell (0-15 Å). A double-peak feature characterizes the IL in
which the anions are Cl-, BF${}_{4}^{-}$ and PF${}_{6}^{-}$, with the two
maxima centered at $\sim 5$ and $\sim 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.NOMAD

###### Acknowledgements.

We gratefully acknowledge computational resources at the HPC Centre of the
University of Strasbourg funded by the Equipex Equip@Meso project (Programme
Investissements d’Avenir) and CPER Alsacalcul/Big Data, and at Grand
Equipement National de Calcul Intensif (GENCI) under allocation DARI
A0100906092 and A0120906092. G.O. and K.F. acknowledge the European Campus
Eucor initiative for Seed Money funding (Project MEDIA). Funding: I.E.A. is
funded by the Idex Contrats Doctoraux – Programme Doctoral International (N.
ESRH0908292D) of the University of Strasbourg. We thank Laurent Douce
(University of Strasbourg) and Carlo Massobrio (University of Strasbourg),
Emanuele Orgiu (INRS, Canada) and Jocelyne Levillain (University of Caen) for
fruitful discussions.

## References

  * (1) W. Chenlu, Y. Wang, Z. Gan, Y. Lu, C. Qian, F. Huo, H. He and S. Zhang, _Chem. Sci._ , 2021, 12, 15503.
  * (2) A. Vioux and B. Coasne, _Adv. Energy Mater._ , 2017, 7, 1700883.
  * (3) T. Welton, _Coord. chem. rev._ , 2004, 248, 2459.
  * (4) J.T. Ye, Y.J. Zhang, R. Akashi, M.S. Bahramy, R. Arita and Y. Iwasa _Science_ , 2012, 338, 1193–1196.
  * (5) S. Jo, D. Costanzo, H. Berger and A. F. Morpurgo, _Nano. Lett._ , 2015, 15, 1197–-1202.
  * (6) S. Jo, D. Costanzo, H. Berger and A. F. Morpurgo, _Nature. Nanotech._ , 2016, 11, 339.
  * (7) C. Jia, M. Famili, M. Carlotti, Y. Liu, P. Wang, I. M. Grace, Z. Feng, Y. Wang, Z. Zhao, M. Ding, X. Xu, C. Wang, S.-J. Lee, Y. Huang, R. C. Chiechi, C. J. Lambert and X. Duan, _Sci. Adv._ , 2018, 4, eaat8237.
  * (8) S. Barbosa Martin, N. Balke, W.-Y., Tsai, C. Santato and M.O. Orlandi, _J. Phys. Chem. Lett._ , 2020, 11, 3257-3262.
  * (9) Y. Pei, Y. Zhang, J. Ma, M. Fan, S. Zhang and J. Wang, _Mater. Today Nano_ , 2022, 17, 2588-8420.
  * (10) D.M. Correia, L.C. Fernandes, P.M. Martins, C. Garcia-Astrain, C.M. Costa, J. Reguera and S. Lanceros-Mendez, _Adv. Funct. Mater._ , 2020, 30, 1909736.
  * (11) L. Song-Lin, K. Tsukagoshi, E. Orgiu and P. Samori, _Chem. Soc. Rev._ , 2016, 45, 118-151.
  * (12) H. Oldamur, F. Malberg, T. Welton, B. Kirchner, _Chem. Soc. Rev._ , 2014, 16, 16880–16890.
  * (13) N. Jose and L. Canongia, _Chem. Soc. Rev._ , 2006, 110, 3330–3335.
  * (14) W. Yanlei, H. He, C. Wang, Y. Lu, K. Dong, F. Huo and S. Zhang, _JACS Au_ , 2022, 2, 543–561.
  * (15) D. Kun, S. Zhang and J. Wang, _Chem. Commun._ , 2021, 52, 6744–6764.
  * (16) D. Bogdan, Q. Yu, N. HC Lewis, W.B. Carpenter, J.M. Bowman and A. Tokmakoff, _Science_ , 2021, 371, 160–164.
  * (17) A.P. Hunt, C.R. Ashworth and R.P. Matthew, _Chem. Soc. Rev._ , 2015, 44, 1257–1288.
  * (18) J.S. Wilkes, _J. Mol. Cat. A: Chem._ , 2004, 214,11-17.
  * (19) P.H. Nguyen, H.N. Duc, K. Hyojung, M.J. Hyung, M.O. Hye and S.J. Mun, _App. Surf. Sci._ , 2022, 596, 153567.
  * (20) M. Michio, S. Shimizu, R. Sotoike, M. Watanabe, Y. Iwasa, Y. Itoh and T. Aida, _J. Am. Chem. Soc._ , 2017, 139, 16072-16075
  * (21) G. Velpula, R. Phillipson, J.X. Lian, D. Cornil, P. Walke, K. Verguts, S. Brems, H. Uji, S. De Gendt, D. Beljonne, R. Lazzaroni, K.S. Mali and S. De Feyter, _ACS Nano_ , 2019, 13, 3512–3521.
  * (22) Y. Yogesh, S.K. Ghosh and S.P. Singh, _ACS Appl. Electr. Mater._ , 2021, 3, 1496-1504
  * (23) P.J. Young, S. Yoo, B. Park, T. Kim, Y.T. Chun, J.M. Kim, K. Kang, S.H. Lee and S.C. Jun, _Nano Research_ , 2021, 14, 2207-2214
  * (24) N. Jose and L. Canongia, _J. Phys. Chem. B_ , 2004, 108, 2038–2047
  * (25) B. Coasne and L. Viau and A. Vioux, _J. Phys. Chem. Lett._ , 2011, 2,1150–1154
  * (26) G. Ori, F. Villemot, L. Viau, A. Vioux and B. Coasne, _Mol. Phys._ , 2014, 112, 1350-1361.
  * (27) G. Ori, C. Massobrio, A. Pradel, M. Ribes and B. Coasne, _Langmuir_ , 2015, 31, 6742-6751.
  * (28) T. Mendez-Morales, J. Carrete, O. Cabeza, L.J. Gallego and L.M. Varela, _J. Phys. Chem. B_ , 2011, 115, 11170-11182.
  * (29) O. Andreussi and N. Marzari, _J. Chem. Phys._ , 2012, 137, 044508.
  * (30) A.A. Niazi, D.B. Rabideau and A.E. Ismail, _J. Phys. Chem. B_ , 2013, 117, 1378-1388.
  * (31) S. Bhattacharjee and S. Khan, _Phys. Chem. Chem. Phys._ , 2020, 22, 8595-8605.
  * (32) S. Agrawal and H.K.Kashyap, _J. Mol. Liq._ , 2018, 261, 337-349.
  * (33) W. Yong-Lei, A. Laaksonen and M.D. Fayer, _J. Phys. Chem. B_ , 2017, 121, 7173-7179.
  * (34) V. Kempter and B. Kirchner, _J. Mol. Structure_ , 2010, 972, 22-34.
  * (35) B. Kirchner, F. Malberg, D.S. Firaha and O. Hollóczki, _J. Phys.: Cond. Matter_ , 2015, 27, 463002.
  * (36) T. Cremer, C. Kolbeck, K.R.J. Lovelock, N. Paape, R. Wolfel, P.S. Schulz, P. Wasserscheid, H. Weber, J. Thar, B. Kirchner, F. Maier and H.‐P. Steinruck, _Chem.-A Eur. J._ , 2010, 16, 9018-9033.
  * (37) P.A. Hunt, B. Kirchner and T. Welton, _Chem.-A Eur. J._ , 2006, 12, 6762-6775.
  * (38) S. Grimme, H. Waldemar and B. Kirchner, _Phys. Chem. Chem. Phys._ , 2012, 14, 4875-4883.
  * (39) M.G. Del Popolo, R.M. Lynden-Bell and J. Kohanoff, _J. Phys. Chem. B_ , 2005, 109, 5895-590.
  * (40) MM. Buhl, A. Chaumont, R. Schurhammer and G. Wipff, _J. Phys. Chem. B_ , 2005, 109, 18591-18599.
  * (41) K. Ishisone, G. Ori and M. Boero, _Phys. Chem. Chem. Phys._ , 2022, 24, 9597-9607.
  * (42) T. Ikeda and M. Boero, _J. Chem. Phys._ , 2008, 143, 194510.
  * (43) S. Gehrke, M. von Domaros, R. Clark, O. Holloczki, M. Brehm, T. Welton, A. Luzar and B. Kirchner, _Phys. Chem. Chem. Phys._ , 2022, 24, 9597-9607.
  * (44) F. Dommert, K. Wendler, R. Berger, L. Delle Site and C. Holm, _ChemPhysChem_ , 2018, 13, 1625-1637.
  * (45) T. Köddermann, D. Reith and R. Ludwig, _ChemPhysChem_ , 2013, 13, 1–8.
  * (46) L. Martínez, R. Andrade, E.G. Birgin and J.M. Martínez, _J. Comput. Chem._ , 2009, 30, 2157–2164.
  * (47) S. Nosé, _Mol. Phys._ , 1984, 52, 255–268.
  * (48) S. Nosé, _J. Chem. Phys._ , 1984, 81, 511–519.
  * (49) W.G. Hoover, _Phys. Rev. A: At., Mol., Opt. Phys._ , 1985, 31, 1695–1697.
  * (50) W. Smith, I.T. Todorov and M. Leslie, _Phys. Rev. A: At., Mol., Opt. Phys._ , 2005, 220, 563–566.
  * (51) W. Kohn and L.J. Sham, _Phys. Rev._ , 1965, 140, A1133–A1138.
  * (52) R. Car and M. Parrinello, _Phys. Rev. Lett._ , 1985, 55, 2471–2474.
  * (53) copyright IBM Corp. (1990-2022), copyright MPI-FKF Stuttgart (1997-2001). www.cpmd.org, v. 4.3, 2022.
  * (54) N. Troullier and J. L. Martins, _Phys. Rev. B: Condens. Matter Mater.Phys._ , 1991, 43, 1993–2006.
  * (55) N. Troullier and J. L. Martins, _Phys. Rev. B: Condens. Matter Mater.Phys._ , 1991, 43, 8861–8869.
  * (56) A. D. Becke, _Phys. Rev. B_ , 1988, 38, 3098.
  * (57) W. Chenlu, Y. Wang, Z. Gan, Y. Lu, C. Qian, F. Huo, H. He, and S. Zhang, _Chem. Sci._ , 2021, 12, 15503.
  * (58) C. Lee, W. Yang and R. G. Parr, _Phys. Rev. B_ , 1988, 37, 785.
  * (59) N. Marzari, A.A. Mostofi, J.R. Yates, I. Souza and D. Vanderbilt, _Rev. Mod. Phys._ , 2012, 84, 1419–1475.
  * (60) P. L. Silvestrelli, _Phys. Rev. Lett._ , 2008, 100, 053002.
  * (61) G. J. Martyna, M. E. Tuckerman and M. L. Klein, _J. Chem. Phys._ , 1992, 97, 2635.
  * (62) P. Söderlind and D.A. Young, _Computation_ , 2018, 6, 13.
  * (63) Nomad Data repository, nomadrepository.eu.