# Modelling proton transfer in [HEIM][TFSI] ionic liquid

Martín Otero-Lema, Kateryna Goloviznina, Luis Miguel Varela, Mathieu Salanne, Hadrian Montes-Campos, Alessandra Serva

► To cite this version:

Martín Otero-Lema, Kateryna Goloviznina, Luis Miguel Varela, Mathieu Salanne, Hadrian Montes-Campos, et al. Modelling proton transfer in [HEIM][TFSI] ionic liquid. Materials today energy, 2025, 53, pp.102018. 10.1016/j.mtener.2025.102018. hal-05343681

HAL Id: hal-05343681

https://hal.science/hal-05343681v1

Submitted on 3 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Modelling proton transfer in [HEIM][TFSI] ionic liquid

Martín Otero-Lema<sup>a,b</sup>, Kateryna Goloviznina<sup>c,d</sup>, Luis Miguel Varela<sup>a,b</sup>, Mathieu Salanne<sup>c,d,e</sup>, Hadrian Montes-Campos<sup>a,b,*</sup>, Alessandra Serva<sup>c,d,** ID</sup>

## ARTICLE INFO

### Keywords:

Proton transfer

Protic ionic liquids

Molecular dynamics

Neural network force field

Polarizable force field

## ABSTRACT

Protic ionic liquids, PILs, are promising materials for energy storage applications, in part due to their ability to decouple proton transport from ion diffusion. In this work, we model the proton transfer mechanism in 1-ethyl-imidazolium bis(trifluoromethanesulfonyl)imide ([HEIM][TFSI]) IL by means of Neural Network Force Field simulations. The latter are combined with classical polarizable molecular dynamics simulations to explore the structure and dynamics of the fully ionized system and Density Functional Theory calculations to estimate the energy barriers for the different proton transfer reactions. Our results show that proton transfer is indeed possible when doping the ionic liquid with an excess of deprotonated cations, but not with an excess of protonated anions. We highlight the importance of the formation of dimers between donor and acceptor species for the reaction to occur, and we identify the main driving factor for the reaction to be the energy cost for reaching a suitable coordination environment and form such dimers, which is higher than that for the transfer reaction.

## 1. Introduction

The ability to store energy is essential for meeting the ever-increasing demands of modern society. To transition from fossil fuels to sustainable energy sources, electrochemical energy storage becomes particularly crucial in ensuring a steady power supply while relying on renewable energy sources such as wind or solar, which are intermittent by nature. Nowadays, the use of electrochemical devices such as lithium-ion batteries [1] (LIBs) and supercapacitors [2] (SCs) is widespread. However, these devices can be hazardous due to the volatility and flammability of the electrolytes needed for their operation [3,4]. Significant efforts have been made to research alternative electrolytes that do not suffer from these problems, while retaining or improving the energy storage capabilities of classical electrolytes.

Ionic Liquids (ILs) have emerged as one of the leading candidates for next-generation electrolytes [5–7]. ILs are composed entirely of cations and anions and can remain liquid at room temperature due to the

difference in size between the ions, their flexibility and the different length of ionic side chains. Together, these factors hinder tight packing, and thus allow thermal vibrations to overcome the Coulomb forces holding them together. ILs exhibit very desirable properties as electrolytes, such as negligible vapor pressure due to the high electrostatic forces between their components, and high thermal stability. These properties make them interesting candidates for novel low-risk ILs or SCs [8–12]. Moreover, the wide array of available cations and anions allows for a high degree of customizability, enabling to fine-tune the properties of the electrolyte. Nevertheless, ILs suffer from high viscosities and low diffusivities, which can hinder their applicability in electrochemical devices, where fast charge transport in the electrolyte is desired [13].

In recent years, much effort has been devoted to the experimental characterization of Protic Ionic Liquid (PIL)-based electrochemical devices [14,15]. Special attention has been paid to studying the interplay between the electrode and the electrolyte, and the phenomenon of

<footnote><sup>a</sup> Grupo de Nanomateriais, Fotónica e Materia Branda, Departamento de Física de Partículas, Universidade de Santiago de Compostela, Campus Vida s/n, Santiago de Compostela, E-15782, Spain</footnote>

<footnote><sup>b</sup> Instituto de Materiais (iMATUS), University of Santiago de Compostela, Campus Vida s/n, Santiago de Compostela, E-15782, Spain</footnote>

<footnote><sup>c</sup> Sorbonne Université, CNRS, Physicochimie des Électrolytes et Nanosystèmes Interenciaux, Paris, F-75005, France</footnote>

<footnote><sup>d</sup> Réseau sur le Stockage Electrochimique de l'Energie (RS2E), FR CNRS 3459, Amiens, 80039, France</footnote>

<footnote><sup>e</sup> Institut Universitaire de France (IUF), Paris, 75231, France</footnote>

<footnote>This article is part of a special issue entitled: AI4EnergyMat published in Materials Today Energy.</footnote>

<footnote>* Corresponding author. Grupo de Nanomateriais, Fotónica e Materia Branda, Departamento de Física de Partículas, Universidade de Santiago de Compostela, Campus Vida s/n, Santiago de Compostela, E-15782, Spain.</footnote>

<footnote>** Corresponding author. Sorbonne Université, CNRS, Physicochimie des Électrolytes et Nanosystèmes Interenciaux, Paris, F-75005, France.</footnote>

<footnote>E-mail addresses: hadrian.montes@usc.es (H. Montes-Campos), alessandra.serva@sorbonne-universite.fr (A. Serva).</footnote>
pseudocapacitance. Materials such as RuO₂, VN and TiN have been found to favour proton transfer between the PIL and the surface [16–18]. Activated carbons have also received attention, since the porous cavities present in these materials can enhance both electrolyte wetting and capacitance [19], while the introduction of functional groups into the carbon matrix was found to increase pseudocapacitance [20].

Theoretical methods have been used to probe the nanostructure of PILs and to complement experimental findings, providing resolution at molecular level. In previous works, the properties of mixtures of PILs with different solvents were studied and compared to those of aprotic ILs, both in bulk and at interfaces [21–23]. The effect of water on these electrolytes has also been extensively studied, since large aggregations of water near the electrode can cause stability issues [24,25]. Recently, A. Balducci et al. have investigated the properties of a novel water-in-PIL electrolyte both experimentally and by means of classical Molecular Dynamics (MD), showing that despite the hygroscopic nature of PILs, it is possible to avoid unwanted reactions at the electrode [26,27].

When it comes to the simulation of PILs, recent developments in simulation techniques have made it possible to correctly describe the phenomena that govern them. On the one hand, the development of polarizable force fields allows to better predict the dynamics, and to improve the description of the hydrogen bond network that permeates them [28–31]. The latter is particularly important, since both integer-charge and scaled-charge non-polarizable force fields are known to underestimate the number of hydrogen bonds in the system [31]. On the other hand, describing proton transfer reactions that take place within some PILs have historically necessitated ab initio methods [32], which are usually very limited in size and time scales that can be investigated. Recently, the rise in popularity of neural networks has lead to the development of neural network force fields (NNFFs) that can extend these scales while retaining the characteristic accuracy of ab initio calculations [33].

In this work, we study the proton transfer mechanism in 1-ethylimidazolium bis(trifluoromethanesulfonyl)imide ([HEIM][TFSI]) IL at different length and time scales. The system is similar to the one previously studied by Hoarfrost and coworkers, who found evidence of proton transfer in mixtures of this PIL with imidazole by means of Nuclear Magnetic Resonance [34] (NMR). This was later confirmed by Mosses et al. using ab initio MD [35]. However, the proton transfer process in the neat IL, as well as its mixtures with protonated or deprotonated variants of the anion and cation, respectively, has received little attention so far. Here, we develop a polarizable force field to describe the different chemical species and to analyze the structure and dynamics of the PIL, as well as to shed light on the possible proton transfers in the system. Then, we characterize the energy barriers of the proton transfer reaction by means of quantum mechanical calculations at the Density Functional Theory (DFT) level. Finally, to go beyond these two approaches, we develop a NNFF to explicitly simulate the proton hopping within the IL medium.

## 2. Methods

### 2.1. Polarizable classical MD simulations

Atomistic MD simulations were carried out in OpenMM [36] software, version 7.6. Two systems containing a total of 600 molecules/ion pairs were considered: a purely ionic, with a 1:1 ratio of [HEIM]⁺ cations and [TFSI]⁻ anions, and a 50 % ionic, formed by a 1:1:1:1 ratio of cations, anions and neutral [EIM] and [HTFSI] species. The corresponding structural formulae are given in Fig. 1.

A polarizable force field was first developed for [HEIM]⁺, [EIM], [TFSI]⁻ and [HTFSI] using the CL&Pol approach [37]. For the charged species, the non polarizable parametrization was taken from the CL&P force field [38] and used as a starting point to build the polarizable models, whereas for the neutral molecules a non polarizable model had

<figure><img src="image_4.png" /><figcaption>Fig. 1. Structural formulae of cation, anion, and neutral species with corresponding labels. Please note that we refer to [HEIM]⁺ also as HC⁺, [EIM] as C, [HTFSI] as HA and [TFSI]⁻ as A⁻, respectively.</figcaption></figure>

to be built before introducing polarization effects (see Supporting Information for all the details).

The simulation procedure was the same for all systems. Initial configurations were built using PACKMOL [39] and ftool [40], and force field files for OpenMM were created with pol_openmm [41]. The systems underwent a 10 ns equilibration in the NpT ensemble, followed by another 5 ns equilibration in the NVT ensemble. Finally, 50 ns production runs in the NVT ensemble were carried out using a timestep of 1 fs. Three independent replicas were simulated for each system, starting from different configurations, in order to adequately sample the phase space. The reported properties were calculated as averages over these replicas. During production runs, positions of all particles were recorded every 0.5 ps, while the electric current vector was recorded at every timestep. For all simulations, temperature was kept at a constant value of 298.15 K by means of a temperature grouped Nose-Hoover thermostat [42,43], whereas in the NpT simulations, pressure was held at 1 atm by means of a Monte Carlo barostat [44,45]. Additional non-equilibrium simulations were performed to compute the viscosity. More details on the simulation setup can be found in the Supporting Information.

### 2.2. DFT calculations

Density functional theory calculations at B3LYP/6-311++G(d,p) level were used to quantify the probability of the proton transfer reactions. Dispersive effects were taken into account using the Grimme model for empirical dispersion with a Becke-Johnson damping [46]. Representative molecular pairs were extracted from the classical MD simulations, and the energy profile was calculated. This was done by gradually increasing the donor-proton distance while allowing the rest of the structure to relax at each step [47]. The effect of the medium was taken into account through a mean field approach (see Supporting Information for more details). Even though this approach does not explicitly consider a condensed phase and assumes only a single pathway, it is useful to verify the feasibility of a reaction.

Once the energy profile along the reaction coordinate was calculated, we followed the approach that Jacobi and coworkers previously used for the 1-methylimidazolium acetate IL to estimate the reaction probability [47]. The normalized Maxwell-Boltzmann kinetic energy distribution for the protonated molecule is given by

$$f(E_k) = \beta^2 E_k e^{-\beta E_k} \quad (1)$$

where $E_k$ is the kinetic energy, and $\beta = (k_B T)^{-1}$, with T absolute temperature and $k_B$ Boltzmann's constant. The reaction probability is then given by the probability that a molecule reaches a kinetic energy higher than the energy barrier

$$p = \int_{\Delta E}^{\infty} f(E_k) dE_k \quad (2)$$

where $\Delta E$ is the height of the energy barrier in the direction of the
reaction. The above equation can be integrated directly to yield the following probability

$$p = (1 + \beta\Delta E) e^{-\beta\Delta E} \quad (3)$$

## 2.3. Neural network force field simulations

The NNFF simulations were carried out using the most recent version of the NeuralIL force field [48] using the ResNet architecture with a committee of 10 neural networks. The training was performed using the same methodology reported in Ref. [48] (VELO optimizer [49] with 51 epochs). The NNFF was developed using configurations extracted from polarizable MD simulations with a reduced set of molecules (7 ion pairs and 5 neutral molecules). DFT calculations using GPAW [50] with the vdW-DF exchange and correlation functional [51] and double zeta polarized basis set were carried out to obtain the energies and forces of 500 configurations extracted at evenly spaced intervals from a 10 ns long polarizable MD simulation in the NpT ensemble, with the same pressure and temperature conditions as previously specified. These data were used to train a first-generation NN which was employed to perform 500 short simulations (250 fs each, and from a different starting configuration). From these new simulations other 500 configurations (the last one of each simulation) were extracted and their energy and forces calculated using DFT. Finally, with a total of 1000 configurations a second generation NN was trained independently from the first one, which was then used for the production simulations. Note that this process does not introduce any bias towards proton transfers. However, it was found that proton transfer was explored by the NN in the second training batch. This can be seen in Fig. 2, where we show the histogram of the reaction coordinate in the two training sets for the NNFF. In the system doped with EIM, where proton transfer is possible, the configurations corresponding to the reaction, where the donor-proton and acceptor-proton distances are equal, are naturally found by the first iteration of the NNFF, and are used to further refine the forcefield.

The production simulations were performed on larger systems using JAX MD [52] in the NVT ensemble for 1.3 ns at 448 K with a Nosé-Hoover thermostat [53]. This higher temperature was chosen in order to accelerate the dynamics of the system, and to be able to observe a significant amount of proton jumps. In fact, NNFF simulations carried at 300 K revealed that no proton diffusion occurs at that temperature for at least 500 ps (see section 3.6.1 of the Supporting Information). The initial configurations were generated using polarizable MD in the NpT ensemble, and consisted of 90 ion pairs and 10 neutral molecules (either EIM or HTFSI). To analyze proton transfer during the simulations, vacant cage correlation functions (CACFs) between hydrogen and acceptor atoms were computed as

$$C(t) = \frac{\langle \theta_i(t) \theta_i(0) \rangle}{\langle \theta_i(0) \theta_i(0) \rangle} \quad (4)$$

where $\theta_i$ takes a value of 1 if the distance between the acceptor atom (i) and all of the hydrogens is greater than a given cutoff $r_0$ and 0 otherwise. The cutoff distances were selected based on coordination distances, since they can be interpreted as the equilibrium distances for a given bond.

# 3. Results and discussion

## 3.1. Polarizable classical MD simulations

We considered two systems, a purely ionic, consisting of a 1:1 ratio of [HEIM]$^+$ cations and [TFSI]$^-$ anions, and a 50% ionic system, formed by a 1:1:1:1 ratio of cations, anions and neutral [EIM] and [HTFSI] species. It should be noted that experimentally the high acidity of the proton in [HTFSI] causes it to fully ionize by transferring its proton to [EIM] [54].

<figure><img src="image_5.png" /><figcaption>Fig. 2. Probability density functions of the proton transfer reaction coordinate in the two training sets for the NNFF of the systems doped with EIM (top) and HTFSI (bottom). The solid lines represent the smoothed surface of the histograms. The reaction coordinate is the difference between the proton-acceptor and proton-donor distances, measured along the donor-acceptor axis. More details on this definition can be found in the Supporting Information.</figcaption></figure>

Because of this, the pure ionic system is the one that better represents the experimental conditions, while the 50% ionic one is used solely to gain insight into the coordination between both charged and neutral species. Therefore, the calculated properties of the 50% ionic system should not be compared with experiments, even though some of the structural features discussed below could be observed experimentally by adding an excess of either [EIM] or [HTFSI] to the pure IL.

### 3.1.1. Force field validation

Since the main advantage of polarizable force fields over their non-polarizable counterparts is the accurate prediction of dynamic properties [31], we validate our model for pure [HEIM] [TFSI] by computing viscosity, electric conductivity and ion diffusion coefficients. Details on the calculations can be found in the Supporting Information. The results are reported in Table 1, together with the calculated density.

A very good agreement is obtained for the density between simulation and experiment, with a relative error not exceeding -0.28 %. For dynamic properties, the agreement is lower, but in all cases the same order of magnitude is obtained. This is an expected outcome of the CL&Pol polarisable force field since it is not fine-tuned to reproduce transport (or other) properties, remaining a general and transferable model for the main families of ionic liquids. The deviations in viscosity range from -2.8 to +6.2 %, the diffusion coefficient of the cation is underestimated by -27 %, and electric conductivity uncertainties vary
<caption>**Table 1**</caption>

Calculated density ($\rho$), viscosity ($\eta$), electric conductivity ($\sigma$) and diffusion coefficients ($D$) of pure [HEIM] [TFSI] at 298 K, compared to the experimental values reported in the ILThermo database [58,59]. The values are averaged between the three replicas, while the uncertainty is estimated from the standard deviation between replicas. Diffusion coefficients and Nernst-Einstein conductivities include the Yeh-Hummer correction [60,61]. Conductivities are calculated using the Green-Kubo (¹) [62], the Einstein-Helfhand (¹) [63,64] and the Nernst-Einstein (³) [65] methods, respectively. Note that the value marked with an asterisk was extrapolated to 298 K using the Arrhenius equation obtained in the 303–353 K range [57]. The numbers in parentheses correspond to the uncertainties in the last decimal places.

<table><thead><tr><th>Property</th><th>Simulation</th><th>Exp</th></tr></thead><tbody><tr><td>ρ / g cm<sup>-3</sup></td><td>1.5685(25)</td><td>1.572(93) [55]<br>1.569(13) [56]</td></tr><tr><td>η / Pa s</td><td>0.0583(16)</td><td>0.0600(24)[55]<br>0.0549(12) [56]<br>0.0576* [57]</td></tr><tr><td>σ / S m<sup>-1</sup></td><td>0.270(18)<sup>†</sup><br>0.279(43)<sup>†</sup><br>0.549(22)<sup>†</sup></td><td>0.368(15) [55]<br>0.244(11) [56]</td></tr><tr><td>D<sub>cat</sub> / 10<sup>-11</sup> m<sup>2</sup> s<sup>-1</sup></td><td>1.84(10)</td><td>2.51(83) [56]</td></tr><tr><td>D<sub>an</sub> / 10<sup>-11</sup> m<sup>2</sup> s<sup>-1</sup></td><td>1.68(10)</td><td>—</td></tr></tbody></table>

from -27 to +14 %. Regarding the latter, three different methods have been used to compute it: the Einstein-Helfhand (EH), the Green-Kubo (GK) and the approximate Nernst-Einstein (NE) methods (see Supporting Information for the details). The comparison between the EH and GK calculations with the approximate NE relation reveals that anion-cation correlations within the liquid contribute significantly to the electric conductivity, lowering its value by approximately a factor of two when compared with the ideal (NE) case. Finally, it is interesting to remark that the obtained GK and EH values are in good agreement with each other, up to statistical uncertainty. When neutral [EIM] and [HTFSI] are added, both viscosity and density decrease, while diffusion coefficients and conductivity increase (see Table S7 of the Supporting Information). This is expected as charged species are substituted with their neutral counterparts, and the system becomes less bound via Coulomb interactions, which should lead to faster dynamics, a well-known effect in ILs [66,67]. In conclusion, the predicted dynamics align well with the experimental data, and the effect of adding neutral molecules match expectations, confirming the force field's suitability for modeling these mixtures.

### 3.1.2. Structural organization of the liquid

In order to study the structure of the system, we calculated several radial distribution functions (RDFs) for both the pure ionic and 50 % ionic systems. The RDFs corresponding to the cation-anion interaction in Fig. 3 (top) and in Fig. S11 of the Supporting Information show that the coordination structure remains mostly unchanged when neutral species are added to the system. In particular, the donor nitrogen atom of [HEIM]<sup>+</sup> has a greater affinity for the oxygen atoms of [TFSI]<sup>-</sup> than for the nitrogen atom. Moreover, as it can be seen in the inset of Fig. 3(top), the coordination between the cation and the O atoms in [TFSI]<sup>-</sup> is one-to-one (i.e the [TFSI]<sup>-</sup> coordinates the cation in a monodentate way), since at the first minimum of the g(r), the coordination number is less than one in both in the pure ionic and 50 % ionic systems.

In Fig. 3 (bottom), the RDFs involving atoms of the neutral molecules are displayed. Concerning [EIM] molecules, their nitrogen atoms, N<sub>3</sub><sup>C</sup>, have strong affinities for the donor nitrogen atoms of both [HEIM]<sup>+</sup> cations and [HTFSI] molecules, hinting at the possible formation of hydrogen bonds between these species. Regarding the coordination between [HTFSI] and [TFSI]<sup>-</sup>, the nitrogen atom in the neutral species, N<sub>3</sub><sup>HA</sup>, shows very different coordinations with the oxygen and nitrogen atoms of the anion. In the case of oxygen atoms, a short range of coordination distances is found, while in the case of the nitrogen atom, a wider range is observed. This suggests that the most probable

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Radial distribution functions (RDFs, g(r)) between characteristic atoms in the pure ionic and 50 % ionic systems. The atom labels are given in Fig. 1, and we refer to [HEIM]<sup>+</sup> as HC<sup>+</sup>, [EIM] as C, [HTFSI] as HA and [TFSI]<sup>-</sup> as A<sup>-</sup>, respectively. Note that O<sup>A<sup>-</sup></sup> refers to the four equivalent oxygen atoms of the anion. The inset in the top panel displays the coordination numbers of anion atoms around the reference cation atom.</figcaption></figure>

coordination is with O atoms, and that the broad distribution in the N<sub>3</sub><sup>HA</sup> - N<sup>A<sup>-</sup></sup> RDF corresponds to said coordinations viewed from the reference N atom, as will be further discussed later.

To further investigate the different coordination environments of hydrogen atoms, combined distribution functions (CDFs) between the donor-acceptor atomic distance (x-axis) and the donor-hydrogen-acceptor angle (y-axis) were calculated from the simulation of the 50 % ionic system. The results are shown in Fig. 4, along with representative MD snapshots of each molecular pair considered. For the cation-anion interaction (Fig. 4a and b), the donor atom of [HEIM]<sup>+</sup> forms a hydrogen bond with an oxygen atom of [TFSI]<sup>-</sup>. However, that hydrogen bond is not formed with the nitrogen atom of the anion, which is the reason for the larger coordination distances observed in the RDF of Fig. 3. These two coordination modes can be analyzed in conjunction with previous NMR spectroscopy results by Hoarfrost and coworkers [34] having reported that no proton conduction occurs in the pure IL. This is consistent with our simulations: a N-H···N hydrogen bond between cations and anions, which could lead to the formation of a pair of neutral molecules and thus contribute to proton conduction, was not observed. Furthermore, the N-H···O hydrogen bond, found in the simulations, does not lead to proton transfer, as that would also induce proton conduction.

In the case of [HTFSI] and [EIM] (Fig. 4c) a hydrogen bond is formed between the two species, in contrast with [TFSI]<sup>-</sup>-[HEIM]<sup>+</sup> pair (Fig. 4a), in which the proton is bound to the cation. This hydrogen bond allows
<figure><img src="image_7.png" /><figcaption>Fig. 4. Combined distribution functions between the donor-acceptor distance (x-axis) and the donor-hydrogen-acceptor angle (y-axis) for different molecules and sites, as well as representative molecular pairs taken from the MD simulation. The colormap reflects the occurrence probability. The color range was chosen based on the maximum occurrence for each configuration, with the exception of subfigure (a), where it was kept similar to that of subfigure (b). The donor-acceptor pairs are: a) $N_3^+$ atom in [HEIM]$^+$ and $N^-$ in [TFSI]$^-$, b) $N_3^+$ atom in [HEIM]$^+$ and $O^-$ atoms in [TFSI]$^-$, c) $N^0$ in [HTFSI] and $N_3^0$ atom in [EIM], d) $N_3^+$ atom in [HEIM]$^+$ and $N_3^0$ atom in [EIM], e) $N^0$ in [HTFSI] and $O^-$ in [TFSI]$^-$, f) $N^0$ in [HTFSI] and $N^-$ in [TFSI]$^-$.</figcaption></figure>

for the ionization of a pair of neutral molecules into a cation-anion pair, as will be discussed later. This is again consistent with experimental findings, which indicate that stoichiometric amounts of [EIM] and [HTFSI] ionize each other [54].

The [HEIM]$^+-$[EIM] pair (Fig. 4d) also shows that they form a hydrogen bond, which should enable proton conduction, as observed experimentally when an excess of deprotonated cations are added to the pure IL [34]. Finally, two different coordination modes can be found between [HTFSI] and [TFSI]$^-$. The first one, in Fig. 4e, shows a shorter coordination distance, and it corresponds to the formation of a hydrogen bond between an oxygen atom in the anion and the nitrogen atom in the neutral molecule. On the other hand, the coordination displayed in Fig. 4f is more interesting, and reveals a proton bridge structure between the two moieties, in a wide range of distances between donor and acceptor. Both these coordinations were found in a study by Munson and coworkers using NMR spectroscopy in conjunction with ab initio calculations [68]. Of the two structures, the proton bridge stands out, as it could enable proton transfer between the two species when an excess of [HTFSI] is added to the pure IL. However, the proton bridge structure is the least common of the two, since the high probability region in Fig. 4f corresponds to the N-N distance for the other dimer. The N-N

coordination distance for the proton bridge structure was found to be around 2.7 Å by inspecting the MD snapshots. Even though the probability of such coordination is not zero, it is very low compared to the N-H···O hydrogen bond (see RDF of Fig. 3).

Finally, it is interesting to consider the potential effect that water can have on the formation of [HEIM]$^+-$[EIM] and [HTFSI]-[TFSI]$^-$ dimers, since water solubility in [HEIM][TFSI] reaches about 7.5 wt % [69]. The presence of water in this liquid has been found to enhance ionic conductivity, as well as diffusion coefficients [70,71], but competition for acceptor sites between water molecules and [HEIM]$^+$/[HTFSI] could hinder the formation of the structures that allow for proton transfer. To investigate this, we carried out a polarizable MD simulation of the 50 % ionic system, to which we added a 10 % molar water content. Analysis of the RDFs of water and donor molecules around acceptor sites (Fig. S28 in the Supporting Information), reveals that, while water molecules interact favourably with acceptor sites, this tendency does not disrupt the structure of the [HEIM]$^+-$[EIM] and [HTFSI]-[TFSI]$^-$ dimers. Thus, small molar fraction of water would not impede the configurations that will lead to proton transfers. Nevertheless, a complete analysis of the influence of water on proton transport would require additional consideration, such as the characterization of the proton transfer
reaction between water molecules and acceptor sites, and will be studied in future works.

## 3.2. DFT proton transfer barriers

As discussed earlier, the existence of hydrogen bonds between different moieties in the system hints at the possibility of proton transfer between those species. To identify good candidates for these reactions, we performed static DFT calculations. The corresponding calculated energy profiles, as well as the probabilities for the forward and backward reactions, are shown in Fig. 5. Panel 5a refers to the proton transfer from [HEIM]$^+$ to [EIM]. Since the reactants and products are identical, there is no difference in the energy minima between the forward and

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Energy profiles for the different proton transfer reactions, as well as the corresponding jumping probabilities calculated from Eq. (3). a: [HEIM]$^+$ + [EIM] $\rightleftharpoons$ [EIM] + [HEIM]$^+$, b: [HTFSI]$_C$ + [TFSI]$_T$ $\rightleftharpoons$ [TFSI]$_C$ + [HTFSI]$_T$, c: [HTFSI]$_C$ + [TFSI]$_C$ $\rightleftharpoons$ [TFSI]$_C$ + [HTFSI]$_C$, d: [HTFSI]$_T$ + [TFSI]$_T$ $\rightleftharpoons$ [TFSI]$_T$ + [HTFSI]$_T$, e: [HEIM]$^+$ + [TFSI]$_C$ $\rightleftharpoons$ [EIM] + [HTFSI]$_C$. The C and T subscripts denote the cis- and trans-isomers of [TFSI]$^-$ and [HTFSI], respectively. The initial and final geometries can be found in Fig. S7 of the Supporting Information.</figcaption></figure>

backward reactions. The reaction probabilities have values around 20 %, i.e. the proton transfer between these molecules is possible. In Fig. 5b–d the energy profiles for the proton transfer between the various isomers of [TFSI]$^-$ and [HTFSI] in the bridged structure (see Fig. S7 of the Supporting Information) are reported. In the case of cis-trans and cis-cis configurations, the proton transfer is possible, albeit with lower probabilities than for [HEIM]$^+$/[EIM]. This, together with the scarcity of the proton bridge structure, could explain why the contribution of hopping to proton conduction is greater when adding [EIM] than when adding [HTFSI] [34]. The trans-trans configuration does not lead to proton conduction, as the mobile proton is captured by an oxygen atom in [TFSI]$^-$, rather than by the nitrogen atom (see Fig. S7), and this process is energetically unfavourable. Finally, it is interesting to note that the two stable configurations are very close in energy for the cis-cis (Fig. 5c), while they are noticeably apart for the cis-trans (Fig. 5b), which can be understood in terms of the different bonding energies of the proton to each of the isomers.

Finally, the energy barrier for the transference of a proton from [HEIM]$^+$ to [TFSI]$^-$ in Fig. 5e shows only one minimum, corresponding to the ionized state, and meaning that the product state is not stable. This is consistent both with the coordination between the two ions, which does not involve a hydrogen bond, and with the coordination between neutral [EIM] and [HTFSI], which does. Indeed, the calculated energy profile suggests that as soon as a hydrogen bond is formed between the neutral molecules, a proton transfer reaction happens and leads to the formation of a cation-anion pair, which is no longer hydrogen bonded, making the backward reaction impossible. Thus, mixing neutral molecules in equal amounts will result in a pure IL where no proton transfers take place.

## 3.3. Neural network force field simulations

In order to explicitly simulate the proton transfer, NNFF simulations were performed for two systems, containing an excess of [EIM] or [HTFSI], respectively. In Fig. 6 (top) the RDF between mobile hydrogen atoms and donor/acceptor nitrogen atoms is represented. A sharp peak at a distance around 1 Å, which corresponds to the hydrogen being covalently bound to the reference atom, can be seen for both systems. A secondary coordination distance around 1.5 Å is observed when adding extra [EIM], which, as will be shown later, corresponds to the proton having jumped to the neutral molecule from the cation. For an excess of [HTFSI], a much broader distribution of distances can be observed, ranging from 2 to 5 Å.

To further analyze this difference, RDFs between donor and acceptor sites are displayed in Fig. 6 (bottom). Here a stark contrast between the two systems is observed. In the case of the addition of [EIM], a clear coordination distance of 3 Å can be observed, which agrees with the results of the polarizable simulations (Fig. 3 (bottom) and Fig. 4d), and which corresponds to a N-H···N hydrogen bond between a [HEIM]$^+$ cation and a [EIM] molecule. On the contrary, when adding an excess of [HTFSI], the coordination distances between anions and [HTFSI] are much broader, with no traces of short coordination distances, again following the same tendencies seen in the classical MD simulations. The fact that the shortest coordination distances start at 3 Å shows that all [TFSI]$^-$-[HTFSI] dimers adopt the configuration shown in Fig. 4e, rather than the bridge structure depicted in Fig. 4f and which exhibits a co-ordination distance of 2.7 Å. Thus, between these species a N-H···O hydrogen bond is formed, rather than a N-H···H bond, and therefore no proton transfer reaction is expected to take place. Fig. 7 also shows the comparison of the RDFs between cations and anions predicted by the NNFF with the polarizable force field. A very good agreement is observed for the N<sub>3</sub><sup>HC</sup> - O<sup>A</sup>- coordination. Regarding the N<sub>3</sub><sup>HC</sup> + N<sup>A</sup>- RDF, the NNFF predicts the existence of a shorter coordination, corresponding to a hydrogen bond between cations and anions, while almost no difference is observed for the second coordination distance. The absence of this hydrogen bond in the polarizable simulations could be
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** RDFs between donor/acceptor nitrogen sites ($N_3^{HC^+}$ and $N_3^S$ atoms) and hydrogen atoms (top), and between donor and acceptor sites (bottom) obtained from the NNFF simulations.</figcaption></figure>

due to the empirical nature of the force field, which does not take into account quantum effects, and the higher temperature of the NNFF simulations, which could allow for the exploration of coordinations not accessible at lower temperatures.

In Fig. 8 we represent the CACFs between hydrogen and donor atoms for the system with an excess of [EIM] and for a variety of cutoff distances, ranging from 1.3 Å to slightly larger distances as 1.5 and 2.0 Å. It is clear that proton transfer occurs in the system, since the decrease of the CACF for long times indicates that some of the hydrogen atoms left the N atom to which they were originally bonded and jumped to a vacant acceptor site. Moreover, two different processes can be clearly appreciated: (i) for short times, the choice of cutoff distance noticeably alters the behaviour of the CACF; (ii) for large times, the three curves collapse into one, independently of $r_0$. The first regime corresponds to the proton hopping back and forth between donor and acceptor, and thus it is sensible to the exact definition of the bound state, while the long time limit can be interpreted as the proton having jumped from the donor to the acceptor, and then for both molecules to have diffused apart. This is to be expected, since once a noticeable number of the molecular pairs have been broken, the specific definition of the bound state will no longer be relevant. However, when analyzing the trajectories for the [HEIM][TFSI] + [HTFSI] system, no proton jumps were found to take place, consistent with the absence of the proton bridge structure. This results in a CACF that takes a constant value of 1 for all choices of $r_0$ (see Fig. S23).

To gain some insight into why no proton diffusion is found in the presence of an excess of [HTFSI], even though the potential energy profiles in Fig. 5 show that it should be possible, energy barriers for the proton transfer reactions were determined from the NNFF simulations.

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Comparison between the radial distribution functions obtained with the polarizable force field and with the NNFF simulations. Note that in the case of the solid red line, temperature effects on the RDF have been taken into account to compare the simulations (see Supporting Information for more details).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Vacant cage auto correlation functions for the hydrogen - acceptor pairs. Here acceptors are taken to be both the $N_3^{HC^+}$ and $N_3^S$ atoms, to account for the indistinguishability of the particles.</figcaption></figure>

We define a reaction coordinate given by $\xi = r_a - r_d$, where $r_a/r_d$ is the distance between the acceptor/donor and the hydrogen atom, projected into the axis that connects the donor and the acceptor (see Fig. S8). Then, the probability that a given donor-acceptor pair has a certain value of this reaction coordinate, $p(\xi)$, is computed from the simulation. Finally, the potential of mean force (PMF), $\Delta E(\xi)$, can be computed as a function of the reaction coordinate as
$$ \Delta E(\xi) = -k_B T \log p(\xi) \qquad (5) $$

The probability distribution was sampled directly from the histograms of the reaction coordinate from the NNFF simulations, without using any enhanced sampling methods. Additional details about the calculation of the PMF can be found in the SI. The results are shown in Fig. 9, along with the energy barriers obtained from DFT calculations, expressed in terms of the new reaction coordinate. In the case of excess [EIM], three energy minima, corresponding to stable configurations, can be seen. The ones at $\xi = \pm 2.75$ Å (labeled $A_{\pm}$) correspond to a [HEIM]$^{+}$-[EIM] pair that are adjacent to one another but where the coordination is not yet optimal for a proton transfer reaction to take place. This state is separated by a large energy barrier of ~4 kcal/mol (corresponding to ~4.5 k$_{B}$T at 448 K) from the most stable configuration, with a value of the reaction coordinate of $\pm 0.5$ Å (labeled $C_{\pm}$). This corresponds to a configuration like the one represented in Fig. 4d, where the proton is shared between two nitrogen atoms, although it is closer to one than to another. Between these two states, an intermediate configuration is found ($B_{\pm}$), which could correspond to rearrangements of the N-H-N angle as the two molecules get closer to one another.

The most stable configuration $C_{\pm}$ has the same values of the reaction coordinate as the energy minima of the energy barriers shown in Fig. 5, and it was chosen as a common origin for both the PMF and the energy barriers computed through DFT. The two calculations are in qualitative agreement, for the location of the minima as well as the width of the proton transfer barrier, although a difference in the height is observed.

<figure><img src="image_12.png" /><figcaption>Fig. 9. Potential of mean force and energy profiles obtained from the NNFF simulations and DFT calculations, respectively, for the proton transfer reaction between [EIM] and [HEIM]$^{+}$ (top) and [TFSI]$^{-}$ and [HTFSI] (bottom). The energy origins for each curve were chosen as that of their respective minima.</figcaption></figure>

This is possibly due to the fact that DFT calculations account for the effect of the environment through a mean field approach, whereas in the NNFF simulations it is done explicitly, and that different exchange-correlation functionals were used.

In summary, the mechanism for proton diffusion for the case of excess [EIM] is the following: Once an adjacent donor-acceptor pair overcomes the energy barrier between the $A_{\pm}$ and $C_{\pm}$ states, they become coordinated and a process begins where the proton hops between the two molecules, transitioning back and forth between the $C_{\pm}$ and $C_{-}$ states. At last, if the large barrier is overcome again, the two molecules diffuse apart, once the proton has been transferred. An illustration of the reaction is shown in Fig. S9 of the Supporting Information, as well as a video of the MD simulation. Finally, an important metric to characterize the proton transfer mechanism is the number of jumps that a hydrogen atom undergoes during a certain time. Following the approach of Minami and Jinnouchi [72], we computed the number of different acceptor sites ($N_3$ atoms) visited by each hydrogen in the system as a function of time, $N(t)$. As shown in Fig. S24 of the Supporting Information, roughly 80% of the hydrogen atoms remain bonded to the same site for the whole simulation, while 20% of them undergo a simple jump from one molecule to another, and 1% of these mobile protons visit three acceptor sites (as shown in the video). This confirms that there is proton diffusion within the system, and not only a simple rattling motion between the $C_{\pm}$ states, and although the kinetic $C_{\pm} \rightarrow A_{\pm}$ barriers shown in Fig. 8 are high, they can still be overcome. Once this happens, the molecular pair breaks apart, and the newly created HEIM$^{+}$ cation forms a new dimer with a different EIM molecule, and transfers the proton to it. The low number of proton transfers in the simulation can be attributed both to the relatively small system sizes and short simulation times. In a larger system, more transfers could happen, even though their contribution to proton diffusion remains small when compared to vehicular transport.

Since the existence of [HEIM]$^{+}$-[EIM] dimers is of critical importance for the proton diffusion phenomenon, we investigated the lifetime of these structures through CACFs (more details on the Supporting Information). The results are shown in Fig. S26 of the Supporting Information, and they reveal that the $C_{\pm}$ and $B_{\pm}$ states are short lived, with characteristic times in the range of hundreds of picoseconds. This is to be expected, since transitions between these two states correspond to a slight rearrangement of the two molecules in the dimer, but it does not entail its dissociation. On the other hand, the $A_{\pm}$ states are found to be more resilient, showing only a very slight decay. This reinforces the importance of the high kinetic barrier between the $A_{\pm}$ and $B_{\pm}$ states, which has to be overcome in order for dimers to be created/destroyed. The CACF shows that this process only happens a few times during the simulations, in good agreement with the low percentage of atoms that visit three different acceptor sites.

Finally, it is possible to investigate the different contributions to proton transport within the system. These will be a diffusive contribution, associated to the vehicular transport of H atoms present in the [HEIM]$^{+}$ cations, as well as a contribution that comes from proton hopping. To distinguish them, we followed an approach that mimics the experimental methodology of Hoarfrost and coworkers [34]. We computed the mean squared displacement (MSD) of H atoms and $N_3$ atoms (following the naming scheme in Fig. 1). The results are shown in Fig. 10, where it can be immediately seen that the motion of H atoms is enhanced when compared to $N_3$ atoms. Moreover, using the well-known Einstein relations (see Supporting Information) it is possible to compute separate diffusion coefficients for protons and nitrogen atoms. We obtain values of $(12.25649 \pm 0.00097) \cdot 10^{-11}$ and $(11.5869 \pm 0.0015) \cdot 10^{-11}$ m$^2$ s$^{-1}$ for proton and nitrogen, respectively. This results in a quotient between the two mechanisms of $D_H/D_N = 1.05779 \pm 0.00016$. This indicates that, while the rate of proton transfer reactions is severely gated by the high energetic cost of dimer formation/separation, their contributions to proton transport within the system is significant.

The case of excess [HTFSI], shown in Fig. 9 (bottom), reveals a
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Mean squared displacement for mobile H atoms and donor/acceptor sites. The dashed lines represent the linear fit carried out to estimate the diffusion coefficient.</figcaption></figure>

similar picture, but with a key difference. Indeed, the PMF displays a single stable state for high reaction coordinates, and a large energy barrier (≥8 kcal/mol, or equivalent to ≥ 9 $k_BT$ at 448 K) is found between this stable state and the configurations explored with DFT. This is due to the fact that, as previously mentioned, no proton bridge structures like the one depicted in Fig. 4f are formed in these simulations, and instead a N-H···O hydrogen bond is formed. The energy barrier in the PMF represents the energy cost associated with a shift from such a configuration to the N-H···N bridge structure, which seems to be very large. This is in agreement with the polarizable MD simulations, where such structures were rare, and explains why no proton jumps are observed even though the barrier for the proton transfer is relatively small (since even if that is the case, the configuration needed for this reaction is not easily reachable).

Therefore, the NNFF simulations, along with the results presented in the previous section, reveal that the nature of proton transport in the studied systems depends upon two key factors: (i) the energy barriers of the proton transfer reactions determine whether or not a proton hop can take place once the involved molecules are coordinated and are forming a hydrogen bond; (ii) reaching a suitable coordination environment is also important, since the corresponding kinetic barriers were found to be higher than those for the proton transfer reaction.

## 4. Conclusions

In this work, the proton transfer mechanism in 1-ethylimidazolium bis(trifluoromethylsulfonyl)-imide IL was studied by means of NNFF MD simulations, accompanied by polarizable MD simulations and quantum mechanical DFT calculations to get prior insights on the system.

Classical MD simulations of both pure and 50 % ionic systems were carried out, after extending a polarizable CL&Pol force field to describe the pure IL, protonated anions and deprotonated cations. They reveal the formation of a strong hydrogen bond between [EIM] and [HTFSI], as neutral species tend to ionize each other to form a cation-anion pair. Moreover, hydrogen bonds between [EIM] and [HEIM]$^+$ were found, as well a N-H···H bond between [HTFSI] and [TFSI]$^-$ (even though less probable), which could enable proton transfer when an excess of either neutral species is added to the pure IL. In addition, it was found that water concentrations up to 10 % molar do not interfere with the structure of the liquid, still allowing for the formation of these hydrogen bonds.

DFT calculations were carried out to identify good candidates for

possible proton transfer reactions. Their analysis confirmed that two neutral species undergo irreversible ionization when coming in contact, while the formation of [EIM]-[HEIM]$^+$ or [TFSI]$^-$-[HTFSI] dimers can lead to a proton transfer reaction. Such a transfer is more likely to happen in the excess of [EIM] with energy barriers of about 1.5 and 3 kcal/mol, respectively. Nevertheless, these absolute values are only approximate, as they account for the dimer solvation environment implicitly and are sensitive to the choice of functional and initial geometries. To explicitly study the proton transfer in mixtures of pure IL with an excess of one of the neutral species, a NNFF was developed. Our simulations reveal that proton transfer indeed takes place within the IL + [EIM] mixture, while no jumps were detected when [HTFSI] is added. This discrepancy was resolved through the analysis of the PMFs. For a proton transfer to take place, the involved molecules must overcome an energy barrier to form a dimer. This is possible in the case of the [HEIM]$^+$-[EIM] pair, for which a barrier of only 4 kcal/mol is required, but not for [HTFSI]-[TFSI]$^-$, where the barrier reaches 9 kcal/mol.

All together, these findings allow to identify the two metrics driving the proton transfer mechanism. On the one hand, the height of the proton transfer barrier between one molecule to another determines if a proton jump is possible once the donor and the acceptor are coordinated and form a dimer. On the other hand, the coordination barrier that the molecules must overcome to form such dimer is also of great importance, since it can prevent otherwise possible proton transfers. The barriers can be accurately captured with NNFF simulations, which achieve the same level of accuracy as DFT results.

Finally, while the developed polarizable force field by its nature cannot take into account proton exchange, it allows for simulating larger sizes and time scales than the NNFF, enabling better exploration of the configuration space and access to dynamic properties. In the future these classical MD simulations could be improved by incorporating proton transfer events, based on DFT reaction probabilities, as it is possible to do with the protex tool [73,74]. This could allow access to rare configurations that are not explored during the NNFF simulations, such as the [HTFSI]-[TFSI]$^-$ proton bridge structure, which will be investigated in upcoming works.

### CRediT authorship contribution statement

**Martín Otero-Lema:** Writing – review & editing, Writing – original draft, Visualization, Investigation, Formal analysis, Data curation.

**Kateryna Goloviznina:** Writing – review & editing, Writing – original draft, Investigation, Formal analysis, Data curation. **Luis Miguel Var-ela:** Writing – review & editing, Supervision, Methodology, Funding acquisition, Conceptualization.

**Mathieu Salanne:** Writing – review & editing, Funding acquisition, Conceptualization. **Hadrian Montes-Campos:** Writing – original draft, Investigation, Formal analysis, Data curation.

**Alessandra Serva:** Writing – review & editing, Writing – original draft, Supervision, Methodology, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

The financial support of the Spanish Ministry of Science and Innovation (PID2021-126148NA-100 funded by MICIU/AEI/10.13039/501100011033/FEDER, UE) is gratefully acknowledged.

This work was funded by the Xunta de Galicia (GRC ED431C 2024/06). M. O. L. thanks the Xunta de Galicia for his "Axudas de apoio á etapa predoutoral" grant (ED481A 2022/236). This work was done within the framework of project HI MOV – "Corredor Tecnológico Transfronterizo de Movilidad con Hidrógeno Renovable", with reference
0160_HI MOV_1_E, co-financed by the European Regional Development fund (ERDF), in the scope of Interreg VI A Spain - Portugal Cooperation Program (POCTEP) 2021-2027. This work is part of the project CNS2023-144785, funded by MICIU/AEI/10.13039/501100011033 and the European Union "NextGenerationEU"/PRTR. H. M. C. thanks the USC for his "Convocatoria de Recualificación do Sistema Universitario Español-Margarita Salas" postdoctoral grant under the "Plan de Recuperación Transformación" program funded by the Spanish Ministry of Universities with European Union's NextGenerationEU funds. The authors acknowledge HPC resources granted by GENCI, France (resources of IDRIS, Grant No. 2023-A0140910463).

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.mtener.2025.102018.

### Data availability

DFT calculations were carried out using Gaussian 16, revision D.01 (https://gaussian.com/) [75], Psi4 version, 1.9.1 (https://psicode.org/) [76] and Gpaw, version 22.1.0 (https://gpaw.readthedocs.io/) [77]. Initial configurations for MD were generated with PACKMOL (https://m3g.github.io/packmol/) [39]. MD simulations were performed using OpenMM (https://openmm.org/) [36] and JaxMD (https://github.com/jax-md/jax-md) [78]. MD trajectories were post-processed using MDAnalysis (https://www.mdanalysis.org/) [79,80], Matplotlib (https://matplotlib.org/) [81], Numpy (https:// numpy.org/) [82], kute (https://gitlab.com/nafomat/kute) [62], Gnuplot (http://www.gnuplot.info/) and TRAVIS (http://www.travis-analyzer.de/) [83,84].

The data and input files that support the findings of this study are available at https://doi.org/10.5281/zenodo.14697374.

## References

[1] M. Armand, J.-M. Tarascon, Building better batteries, Nature 451 (2008) 652–657, 10.1038/451652a.

[2] P. Simon, Y. Gogotsi, Perspectives for electrochemical capacitors and related devices, Nat. Mater. 19 (2020) 1151–1163, https://doi.org/10.1038/s41563-020-0747-z.

[3] G. Walden, J. Stepan, C. Mikolajczak, Safety considerations when designing portable electronics with electric double-layer capacitors (supercapacitors), in: 2011 IEEE Symposium on Product Compliance Engineering Proceedings, IEEE, 2011, pp. 1–5, https://doi.org/10.1109/PSES.2011.6088259.

[4] Q. Wang, L. Jiang, Y. Yu, J. Sun, Progress of enhancing the safety of lithium ion battery from the electrolyte aspect, Nano Energy 55 (2019) 93–114, https://doi.org/10.1016/j.jaoen.2018.10.035.

[5] A. Lewandowski, A. Swiderska-Moczek, Ionic liquids as electrolytes for li-ion batteries—an overview of electrochemical studies, J. Power Sources 194 (2) (2009) 601–609, https://doi.org/10.1016/j.jpowsour.2009.06.089.

[6] M. Salanne, Ionic liquids for supercapacitor applications, Top. Curr. Chem. 375 (2017) 63, https://doi.org/10.1007/s41061-017-0150-7.

[7] W. Zhou, M. Zhang, X. Kong, W. Huang, Q. Zhang, Recent advance in ionic-liquid-based electrolytes for rechargeable metal-ion batteries, Adv. Sci. 8 (13) (2021) 2004490, https://doi.org/10.1002/advs.202004490.

[8] M. Galiński, A. Lewandowski, I. Stępniaj, Ionic liquids as electrolytes, Electrochim. Acta 51 (26) (2006) 5567–5580, https://doi.org/10.1016/j.electacta.2006.03.016.

[9] D.S. Silvester, R.G. Compton, Electrochemistry in room temperature ionic liquids: a review and some possible applications, Z. Phys. Chem. 220 (10) (2006) 1247–1274, https://doi.org/10.1524/zpch.2006.220.10.1247.

[10] H. Liu, Y. Liu, J. Li, Ionic liquids in surface electrochemistry, Phys. Chem. Chem. Phys. 12 (8) (2010) 1685–1697, https://doi.org/10.1039/B921469K.

[11] K. Ghandi, A review of ionic liquids, their limits and applications, Green Sustain. Chem. 4 (2014) 44–53, https://doi.org/10.4236/gsc.2014.41008.

[12] E. Jónsson, Ionic liquids as electrolytes for energy storage applications—a modelling perspective, Energy Storage Mater. 25 (2020) 827–835, https://doi.org/10.1016/j.ensm.2019.08.030.

[13] Z. Lin, E. Goikolea, A. Balducci, K. Naoi, P.-L. Taberna, M. Salanne, G. Yushin, P. Simon, Materials for supercapacitors: when li-ion battery power is not enough, Mater. Today 21 (2018) 419–436, https://doi.org/10.1016/j.mattod.2018.01.035.

[14] C. Zhao, G. Burrell, A.A. Torriero, F. Separovic, N.F. Dunlop, D.R. MacFarlane, A. M. Bond, Electrochemistry of room temperature protic ionic liquids, J. Phys. Chem. B 112 (23) (2008) 6923–6936, https://doi.org/10.1021/jp711804j.

[15] S. Menne, J. Pires, M. Anouti, A. Balducci, Protic ionic liquids as electrolytes for lithium-ion batteries, Electrochem. Commun. 31 (2013) 39–41, https://doi.org/10.1016/j.elecom.2013.02.026.

[16] S. Ketabi, B. Decker, K. Lian, Proton conducting ionic liquid electrolytes for liquid and solid-state electrochemical pseudocapacitors, Solid State Ionics 298 (2016) 73–79, https://doi.org/10.1016/j.ssi.2016.11.010.

[17] A. Djire, J.Y. Ishimwe, S. Choi, L.T. Thompson, Enhanced performance for early transition metal nitrides via pseudocapacitance in protic ionic liquid electrolytes, Electrochem. Commun. 77 (2017) 19–23, https://doi.org/10.1016/j.elecom.2017.02.001.

[18] J.S. Seenath, D. Pech, D. Rochefort, Investigation of protic ionic liquid electrolytes for porous ruo2 micro-supercapacitors, J. Power Sources 548 (2022) 232040, https://doi.org/10.1016/j.jpowsour.2022.232040.

[19] M. Anouti, E. Couadou, L. Timperman, H. Galiano, Protic ionic liquid as electrolyte for high-densities electrochemical double layer capacitors with activated carbon electrode material, Electrochim. Acta 64 (2012) 110–117, https://doi.org/10.1016/j.electacta.2011.12.120.

[20] R. Mysyk, E. Raymundo-Piñero, M. Anouti, D. Lemordant, F. Béguin, Pseudo-capacitance of nanoporous carbons in pyrrolidinium-based protic ionic liquids, Electrochem. Commun. 12 (3) (2010) 414–417, https://doi.org/10.1016/j.elecom.2010.01.007.

[21] P. Martínez-Crespo, M. Otero-Lema, O. Cabeza, H. Montes-Campos, L.M. Varela, Structure, dynamics and ionic conductivities of ternary ionic liquid/lithium salt/dmso mixtures, J. Mol. Liq. 359 (2022) 119188, https://doi.org/10.1016/j.molliq.2022.119188.

[22] M. Otero-Lema, P. Martínez-Crespo, T. Méndez-Morales, H. Montes-Campos, L. M. Varela, Interfacial structure of protic and aprotic ionic liquid-dmso-li salt mixtures near charged and neutral electrodes: a molecular dynamics study, J. Mol. Liq. 386 (2023) 122492, https://doi.org/10.1016/j.molliq.2023.122492.

[23] R. Lois-Cuns, M. Otero-Lema, A. Rivera-Pousa, P. Vallet, J.J. Parajó, O. Cabeza, H. Montes-Campos, T. Méndez-Morales, L.M. Varela, Mixtures of ethylammonium nitrate and ethylene carbonate: bulk and interfacial analysis, J. Mol. Liq. 385 (2023) 122361, https://doi.org/10.1016/j.molliq.2023.122361.

[24] K. Xu, Electrolytes and interphases in li-ion batteries and beyond, Chem. Rev. 114 (23) (2014) 11503–11618, https://doi.org/10.1021/cr500003w.

[25] B. Yang, A. Hu, T. Li, K. Li, Y. Li, J. Jiang, Z. Xiao, Z.W. Seh, J. Long, Eliminating water hazards and regulating electrode-electrolyte interfaces by multifunctional sacrificial electrolyte additives for long-life lithium metal batteries, Energy Storage Mater. 70 (2024) 103512, https://doi.org/10.1016/j.ensm.2024.103512.

[26] T. Stettner, S. Gehrke, P. Ray, B. Kirchner, A. Balducci, Water in protic ionic liquids: proper ties and use of a new class of electrolytes for energy-storage devices, ChemSusChem 12 (16) (2019) 3827–3836, https://doi.org/10.1002/cssc.201901283.

[27] L. Dick, T. Stettner, Y. Liu, S. Liu, B. Kirchner, A. Balducci, Hygroscopic protic ionic liquids as electrolytes for electric double layer capacitors, Energy Storage Mater. 53 (2022) 744–753, https://doi.org/10.1016/j.ensm.2022.09.025.

[28] K. Goloviznina, Z. Gong, A.A.H. Padua, The CL&P polarizable force field for the simulation of ionic liquids and eutectic solvents, Wiley Interdiscip. Rev. Comput. Mol. Sci. 12 (3) (2022) e1572, https://doi.org/10.1002/wcm5.1572.

[29] F. Joerg, C. Schröder, Polarizable molecular dynamics simulations on the conductivity of pure 1-methylimidazolium acetate systems, Phys. Chem. Chem. Phys. 24 (25) (2022) 15245–15254, https://doi.org/10.1039/D2CP01501C.

[30] F. Joerg, J. Sutter, L. van Dam, K. Kanellopoulos, J. Hunger, C. Schröder, Comparative analysis of dielectric spectra in protic ionic liquids: experimental findings and computational molecular decomposition, J. Mol. Liq. 396 (2024) 123834, https://doi.org/10.1016/j.molliq.2023.123834.

[31] D. Bedrov, J.-P. Piquemal, O. Borodin, A.D. MacKerell, B. Roux, C. Schröder, Molecular dynamics simulations of ionic liquids and electrolytes using polarizable force fields, Chem. Rev. 119 (13) (2019) 7940–7995, https://doi.org/10.1021/acs.chemrev.8b00763.

[32] A. Le Donne, E. Bodo, Isomerization patterns and proton transfer in ionic liquids constituents as probed by ab-initio computation, J. Mol. Liq. 249 (2018) 1075–1082, https://doi.org/10.1016/j.molliq.2017.11.152.

[33] H. Montes-Campos, J. Carrete, S. Bichelmaier, L.M. Varela, G.K. Madsen, A differentiable neural-network force field for ionic liquids, J. Chem. Inf. Model. 62 (1) (2021) 88–101, https://doi.org/10.1021/acs.jcim.1c01380.

[34] M.L. Hoarfrost, M. Tyagi, R.A. Segalman, J.A. Reimer, Proton hopping and long-range transport in the protic ionic liquid [im][tfsi], probed by pulsed-field gradient nmr and quasi-elastic neutron scattering, J. Phys. Chem. B 116 (28) (2012) 8201–8209, https://doi.org/10.1021/jp3044237.

[35] A.A. Moses, C. Arntsen, Ab initio molecular dynamics study of proton transport in imidazolium-based ionic liquids with added imidazole, Phys. Chem. Chem. Phys. 25 (3) (2023) 2142–2152, https://doi.org/10.1039/D2CP03262G.

[36] P. Eastman, J. Swails, J.D. Chodera, R.T. McGibbon, Y. Zhao, K.A. Beauchamp, L.-P. Wang, A.C. Simmonett, M.P. Harrigan, C.D. Stern, et al., Openmm 7: rapid development of high performance algorithms for molecular dynamics, PLoS Comput. Biol. 13 (7) (2017) e1005659, https://doi.org/10.1371/journal.pcbi.1005659.

[37] K. Goloviznina, J.N. Canongia Lopes, M. Costa Gomes, A.A.H. Padua, Transferable, polarizable force field for ionic liquids, J. Chem. Theor. Comput. 15 (11) (2019) 5858–5871, https://doi.org/10.1021/acs.jctc.9b00689.

[38] J.N. Canongia Lopes, A. Pádua, CL&P: a generic and systematic force field for ionic liquids modeling, Theor. Chem. Acc. 131 (2012) 1129, https://doi.org/10.1007/s00214-012-1129-7.
[39] L. Martínez, R. Andrade, E.G. Birgin, J.M. Martínez, PACKMOL: a package for building initial configurations for molecular dynamics simulations, J. Comput. Chem. 30 (13) (2009) 2157–2164, https://doi.org/10.1002/jcc.21224.

[40] A.A.H. Padua. https://github.com/paduagroup/fftool, 2024.

[41] A.A.H. Padua. https://github.com/paduagroup/pol_openmm, 2024.

[42] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, A smooth particle mesh ewald method, J. Chem. Phys. 103 (19) (1995) 8577–8593, https://doi.org/10.1063/1.470117.

[43] Z. Gong, A.A. Padua, Effect of side chain modifications in imidazolium ionic liquids on the properties of the electrical double layer at a molybdenum disulfide electrode, J. Chem. Phys. 154 (8) (2021) 084504, https://doi.org/10.1063/5.0040172.

[44] K.-H. Chow, D.M. Ferguson, Isothermal-isobaric molecular dynamics simulations with Monte Carlo volume sampling, Comput. Phys. Commun. 91 (1–3) (1995) 283–289, https://doi.org/10.1016/0010-4655(95)00059-0.

[45] J. Åqvist, P. Wennerström, M. Nervall, S. Bjelic, B.O. Brandsdal, Molecular dynamics simulations of water and biomolecules with a Monte Carlo constant pressure algorithm, Chem. Phys. Lett. 384 (4–6) (2004) 288–294, https://doi.org/10.1016/j.cplett.2003.12.039.

[46] S. Grimme, S. Ehrlich, L. Goerigk, Effect of the damping function in dispersion corrected density functional theory, J. Comput. Chem. 32 (7) (2011) 1456–1465, https://doi.org/10.1002/jcc.21759.

[47] R. Jacobi, F. Joerg, O. Steinhauser, C. Schröder, Emulating proton transfer reactions in the pseudo-protic ionic liquid 1-methylimidazolium acetate, Phys. Chem. Chem. Phys. 24 (16) (2022) 9277–9285, https://doi.org/10.1039/D2CP00643J.

[48] J. Carrete, H. Montes-Campos, R. Wanzenböck, E. Heid, G.K. Madsen, Deep ensembles vs committees for uncertainty estimation in neural-network force fields: comparison and application to active learning, J. Chem. Phys. 158 (20) (2023) 204801, https://doi.org/10.1063/5.0146905.

[49] L. Metz, J. Harrison, C.D. Freeman, A. Merchant, L. Beyer, J. Bradbury, N. Agrawal, B. Poole, I. Mordatch, A. Roberts, et al., Velo: training versatile learned optimizers by scaling up, arXiv preprint arXiv:2211.09760 (2022), https://doi.org/10.48550/arXiv.2211.09760.

[50] J. Enkovaara, C. Rostgaard, J.J. Mortensen, J. Chen, M. Dulak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. Hansen, et al., Electronic structure calculations with gpaw: a real-space implementation of the projector augmented-wave method, J. Phys. Condens. Matter 22 (25) (2010) 253202, https://doi.org/10.1088/0953-8984/22/25/253202.

[51] V.R. Cooper, Van der waals density functional: An appropriate exchange functional, Phys. Rev. B Condens. Matter 81 (16) (2010) 161104, https://doi.org/10.1103/PhysRevB.81.161104.

[52] S.S. Schoenholz, E.D. Cubuk, Jax md: end-to-end differentiable, hardware accelerated. Molecular Dynamics in Pure python, 2019.

[53] D.J. Evans, B.L. Holian, The nose-hoover thermostat, J. Chem. Phys. 83 (8) (1985) 4069–4074, https://doi.org/10.1063/1.449071.

[54] T. Yasuda, M. Watanabe, Protic ionic liquids: fuel cell applications, MRS Bull. 38 (7) (2013) 560–566, https://doi.org/10.1557/mrs.2013.153.

[55] T. Makino, M. Kanakubo, Y. Masuda, T. Umecky, A. Suzuki, Co2 absorption properties, densities, viscosities, and electrical conductivities of ethylimidazolium and 1-ethyl-3-methylimidazolium ionic liquids, Fluid Phase Equilib. 362 (2014) 300–306, https://doi.org/10.1016/j.fluid.2013.10.031.

[56] I. Abdurrokhman, K. Elamin, O. Danyliv, M. Hasani, J. Swenson, A. Martinelli, Protic ionic liquids based on the alkyl-imidazolium cation: effect of the alkyl chain length on structure and dynamics, J. Phys. Chem. B 123 (18) (2019) 4044–4054, https://doi.org/10.1021/acs.jpcb.9b01274.

[57] T. Vogl, P. Goodrich, J. Jacquemin, S. Passerini, A. Balducci, The influence of cation structure on the chemical-physical properties of protic ionic liquids, J. Phys. Chem. C 120 (16) (2016) 8525–8533, https://doi.org/10.1021/acs.jpcc.b061945.

[58] Q. Dong, A.F. Kazakov, C.D. Muzny, R.D. Chirico, J.A. Widegren, V. Diky, J. W. Magee, K.N. Marsh, M.D. Frenkel, Ionic Liquids Database (Ilthermo), 2006.

[59] Q. Dong, C.D. Muzny, A. Kazakov, V. Diky, J.W. Magee, J.A. Widegren, R. D. Chirico, K.N. Marsh, M. Frenkel, Ilthermo: a free-access web database for thermodynamic properties of ionic liquids, J. Chem. Eng. Data 52 (4) (2007) 1151–1159, https://doi.org/10.1021/je700171f.

[60] I.-C. Yeh, G. Hummer, System-size dependence of diffusion coefficients and viscosities from molecular dynamics simulations with periodic boundary conditions, J. Phys. Chem. B 108 (40) (2004) 15873–15879, https://doi.org/10.1021/jp0477147.

[61] E.J. Maginn, R.A. Messerly, D.J. Carlson, D.R. Roe, J.R. Elliot, Best practices for computing transport properties 1. self-diffusivity and viscosity from equilibrium molecular dynamics [article v1.0], J. Comp. Mol. Sci. 1 (1) (2019) 6324, https://doi.org/10.33011/livecoms.1.1.6324, 6324.

[62] M. Otero-Lema, R. Lois-Cuns, M.A. Boado, H. Montes-Campos, T. Méndez-Morales, L.M. Varela. https://gitlab.com/nafomat/kute, 2024.

[63] C. Schröder, M. Haberler, O. Steinhauser, On the computation and contribution of conductivity in molecular ionic liquids, J. Chem. Phys. 128 (13) (2008) 134501, https://doi.org/10.1063/1.2868752.

[64] C. Schröder, O. Steinhauser, On the dielectric conductivity of molecular ionic liquids, J. Chem. Phys. 131 (11) (2009) 114504, https://doi.org/10.1063/1.3220069.

[65] H.K. Kashyap, H.V. Annapureddy, F.O. Raineri, C.J. Margulis, How is charge transport different in ionic liquids and electrolyte solutions? J. Phys. Chem. B 115 (45) (2011) 13212–13221, https://doi.org/10.1021/jp204182c.

[66] C.A. Brooks, A.P. Doherty, Concentration-dependent diffusion in room temperature ionic liquids: a microelectrode study, Electrochem. Commun. 6 (8) (2004) 867–871, https://doi.org/10.1016/j.elecom.2004.06.012.

[67] J. Cui, T. Kobayashi, R.L. Sacci, R.A. Matsumoto, P.T. Cummings, M. Pruski, Diffusivity and structure of room temperature ionic liquid in various organic solvents, J. Phys. Chem. B 124 (44) (2020) 9931–9937, https://doi.org/10.1021/acs.jpcb.0c07582.

[68] K.T. Munson, J. Vergara, L. Yu, T.D. Vaden, Characterization of the bridged proton structure in htfsi acid ionic liquid solutions, J. Phys. Chem. B 119 (20) (2015) 6304–6310, https://doi.org/10.1021/acs.jpcb.5b02715.

[69] M.A.R. Martins, C.M.S.S. Neves, K.A. Kurnia, L.M.N.B.F. Santos, M.G. Freire, S. P. Pinho, J.A.P. Coutinho, Analysis of the isomerism effect on the mutual solubilities of bis(trifluoromethylsulfonyl)imide-based ionic liquids with water, Fluid Phase Equilib. 381 (2014) 28–35, https://doi.org/10.1016/j.fluid.2014.08.007.

[70] N. Yaghini, V. Gómez-González, L.M. Varela, A. Martinelli, Structural origin of proton mobility in a protic ionic liquid/imidazolium mixture: insights from computational and experimental results, Phys. Chem. Chem. Phys. 18 (33) (2016) 23195–23206, https://doi.org/10.1039/C6CP03058K.

[71] M. Hasani, L.M. Varela, A. Martinelli, Short-range order and transport properties in mixtures of the protic ionic liquid [c2him] [tfsi] with water or imidazole, J. Phys. Chem. B 124 (9) (2020) 1767–1777, https://doi.org/10.1021/acs.jpcb.9b10454.

[72] S. Minami, R. Jinnouchi, Accelerating anhydrous proton conduction via anion rotation and hydrogen bond recombination: a machine-learning molecular dynamics, J. Mater. Chem. A 11 (30) (2023) 16104–16114, https://doi.org/10.1039/D3TA03164K.

[73] F. Joerg, M. Wieder, C. Schröder, Protex—a python utility for proton exchange in molecular dynamics simulations, Front. Chem. 11 (2023) 1140896, https://doi.org/10.3389/fchem.2023.1140896.

[74] M. Gödeny, F. Joerg, M.P.-P. Kovar, C. Schröder, Updates to protex for simulating proton transfers in an ionic liquid, J. Phys. Chem. B 128 (14) (2024) 3416–3426, https://doi.org/10.1021/acs.jpcb.3c07356.

[75] M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmanvi, B. Barone, G.A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B.G. Janesko, R. Gomperts, B. Mennucci, H.P. Hratchian, J. V. Ortiz, A.F. Iz- maylov, J.L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J.A. Montgomery Jr., J.E. Peralta, F. Ogliaro, M. J. Bearpark, J.J. Heyd, E.N. Brothers, K.N. Kudin, V.N. Staroverov, T.A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A.P. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, J.M. Millam, M. Klene, C. Adamo, R. Cammi, J.W. Ochterski, R.L. Martin, K. Morokuma, O. Farkas, J.B. Foresman, D.J. Fox, Gaus- Sian16 Revision B.01, gaussian Inc., Wallingford CT, 2016. URL, https://gaussian.com/gaussian16/.

[76] D.G. Smith, L.A. Burns, A.C. Simmonett, R.M. Parrish, M.C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, et al., Psi4 1.4: open-source software for high-throughput quantum chemistry, J. Chem. Phys. 152 (18) (2020) 184108, https://doi.org/10.1063/5.0006002.

[77] J.J. Mortensen, A.H. Larsen, M. Kuisma, A.V. Ivanov, A. Taghizadeh, A. Peterson, A. Haldar, A.O. Dohn, C. Schäfer, E.O. Jónsson, et al., Gpaw: an open python package for electronic structure calculations, J. Chem. Phys. 160 (9) (2024) 092503, https://doi.org/10.1063/5.0182685.

[78] S.S. Schoenholz, E.D. Cubuk, Jax m.d. a framework for differentiable physics, in: Advances in Neural Information Processing Systems, vol. 33, Curran Associates, Inc., 2020. URL, https://papers.nips.cc/paper/2020/file/83d3d4b6c9579515e1679aca8cbc3033-Paper.pdf.

[79] N. Michaud-Agrawal, E.J. Denning, T.B. Woolf, O. Beckstein, Mdanalysis: a toolkit for the analysis of molecular dynamics simulations, J. Comput. Chem. 32 (10) (2011) 2319–2327, https://doi.org/10.1021/jc217788n.

[80] R.J. Gowers, M. Linke, J. Barnoud, T.J.E. Reddy, M.N. Melo, S.L. Seyler, J. Do-manski, D.L. Dotson, S. Buchoux, I.M. Kenney, et al., Mdanalysis: a python package for the rapid analysis of molecular dynamics simulations, Tech. rep., Los Alamos Na- tional Laboratory (LANL), Los Alamos, NM, 2019, https://doi.org/10.25080/Majora-629e541a-00e.

[81] J.D. Hunter, Matplotlib: a 2d graphics environment, Comput. Sci. Eng. 9 (3) (2007) 90–95, https://doi.org/10.1109/MCSE.2007.55.

[82] C.R. Harris, K.J. Millman, S.J. Van Der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N.J. Smith, et al., Array programming with numpy, Nature 585 (7825) (2020) 357–362, https://doi.org/10.1038/s41586-020-2649-2.

[83] M. Brehm, B. Kirchner, Travis - a free analyzer and visualizer for Monte Carlo and molecular dynamics trajectories, J. Chem. Inf. Model. 51 (8) (2011) 2007–2023, https://doi.org/10.1021/ci200217w.

[84] M. Brehm, M. Thomas, S. Gehrke, B. Kirchner, Travis—a free analyzer for trajectories from molecular simulation, J. Chem. Phys. 152 (16) (2020) 164105, https://doi.org/10.1063/5.0005078.