Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Molecular-level heterogeneity in deep eutectic electrolytes

Mirna Alhanash a, Carolina Cruz a and Patrik Johansson *ab
aDepartment of Physics, Chalmers University of Technology, 412 96 Göteborg, Sweden. E-mail: patrik.johansson@chalmers.se
bAlistore-ERI, CNRS FR 3104, 15 Rue Baudelocque, 80039 Amiens, France

Received 15th May 2025 , Accepted 26th August 2025

First published on 5th September 2025


Abstract

Deep eutectic electrolytes (DEEs) have in recent years gained momentum as prospective compounds for next-generation (lithium) batteries. Yet, the connection between molecular-level properties and macroscopic performance of DEEs is poorly understood. We have here, by molecular dynamics simulations, studied in detail a few simple DEEs created using the hydrogen bond (HB) donor N-methyl-acetamide (NMA) and one of three different common lithium salts: LiBF4, LiDFOB, and LiBOB, all in a 1[thin space (1/6-em)]:[thin space (1/6-em)]4 molar ratio composition. By altered anion size and symmetry, we observe varying molecular-level heterogeneity (MLH), reflecting differences in local structure and coordination as well as dynamic disorder, and especially the heterogeneity of the HB network. The HB network becomes more localized the more asymmetric and larger the anion is, which also affects the ion-pairing and aggregation of the DEEs. From a dynamic point of view, DEEs with larger MLH show slower ion self-diffusion, due to steric hindrances caused by the larger anions and the localized HB network. Overall, both the MLH and the HB network, inherent properties of DEEs, must be properly balanced to allow for adequate ion transport.


Introduction

The growing need for renewable energy integration and transport electrification requires efficient energy storage, i.e., battery technologies.1 The electrolyte plays a crucial role for any battery technology, and while liquid electrolytes have been extensively researched for lithium-ion batteries (LIBs), they are limited by e.g., the thermal stabilities of the solvents and salts employed.2,3 Deep eutectic electrolytes (DEEs) have emerged as a promising alternative to conventional liquid electrolytes and have been used both for LIBs and lithium metal batteries (LMB).2,4–8 DEEs-in this context-are eutectic mixtures of pure compounds, often based on the sole addition of a hydrogen bond donor (HBD) to a lithium salt.2,9,10 DEEs offer several advantages, such as ease of synthesis, tunability, wide electrochemical stability windows, and low vapor pressures.2,11

DEEs have, for example, been applied to create LMBs using lithium bis(trifluoromethanesulfonyl) imide (LiTFSI) in combination with 2,2′-dipyridyl disulfide (DpyDS) as HBD, showing a strong ion–dipole interaction promoting enhanced ion conductivity.12 Additionally, Hu et al.13 developed a DEE based on LiTFSI and acetamide, claiming exceptional physicochemical properties with an ionic conductivity exceeding 1 mS cm−1. We here, however, focus on three different lithium salts, whereof lithium difluoro(oxalato)borate (LiDFOB) based electrolytes render improved ionic conductivities as compared to lithium tetrafluoroborate (LiBF4) based ones, and in addition, a better-working low-temperature working range for LIBs.14–16 Our third salt of choice, lithium bis(oxalato)borate (LiBOB), has gained interest mainly due to its fluorine-free design, lower cost, and higher thermal stability as compared to the for LIBs commonly used LiPF6 salt.17,18

However, the main reason as to why we here have chosen these three salts is that they represent a variation in anion symmetry and size (Table 1). The BF4 anion, the highest in symmetry, is pseudospherical, while formally belonging to the Td point group, and is also the smallest. Moving to bulkier and lower in symmetry anions, BOB is elongated and has rings rich in oxygen sites, in stark contrast to BF4 with no such sites at all. DFOB naturally presents in-between characteristics.

Table 1 Anion structure, symmetry, and size
Anion Structure Symmetry # Symmetry elements Sizea
a Length of longest axis = maximum atomic distance.
BF4 image file: d5cp01828e-u1.tif T d 20 ∼2.3 Å
DFOB image file: d5cp01828e-u2.tif D 2h 8 ∼4.3 Å
BOB image file: d5cp01828e-u3.tif D 2d 7 ∼6.6 Å


Overall, these differences in symmetry, size, and chemical composition render different charge distributions, as visualized by their electrostatic potentials (ESPs) (Fig. 1). This strategy thus allows us to systematically monitor effects on DEE properties.


image file: d5cp01828e-f1.tif
Fig. 1 ESPs for: (a) BF4, (b) DFOB, (c) BOB, (d) NMA.

As the HBD, we have chosen N-methyl acetamide (NMA), for which, e.g. the work by Boisset et al.,5 shows that it acts as a complexing agent for both anions and cations and can thus potentially disrupt the cation–anion interactions. This functionality arises from its two polar groups—the NHR group, which coordinates with anions, and the CO group, which coordinates with cations–visible by its ESP (Fig. 1).

Bang and Yeu19 also highlighted the influence of the composition, i.e., the Li-salt:NMA molar ratio, on the ionic and electronic conductivity. Their DEEs exhibited increased ionic conductivity with decreased salt concentration when moving from molar ratios 1[thin space (1/6-em)]:[thin space (1/6-em)]1 to 1[thin space (1/6-em)]:[thin space (1/6-em)]5. While there is no fixed ratio that guarantees the eutectic point, the optimal molar ratio for ionic conductivity was between 1[thin space (1/6-em)]:[thin space (1/6-em)]4 to 1[thin space (1/6-em)]:[thin space (1/6-em)]5.

The more general class of deep eutectic solvents (DESs), which form the basis of DEEs, have been extensively studied for a wide range of applications,20,21 including as battery electrolytes, and their fundamental characteristics and behaviors carry more or less directly over into DEEs.4–8 DESs have recently been reported to have a specific molecular ordering and a radially layered structure,22,23 which is suggested to arise due to a combination of interactions, such as Coulombic, hydrogen bonding (HB), and van der Waals interactions, etc.,22–24 and has also been described in terms of spatial heterogeneity on the molecular length scale.20 Several of these insights have been obtained using MD simulations, highlighting the role of HB networks and structural organization in DESs systems.20,22–25 One study described what was referred to as micro-heterogeneity (MH) as the separation of polar and nonpolar molecular parts, revealing that the HB network of the polar domains could be indirectly affected by the MH.25 A complex HB network seems to be a main characteristic of DESs,20,23–29 and we hypothesize here that this is the case also for DEEs. However, DEEs remain less studied from a molecular simulation perspective, and the specific role of anion symmetry and size in shaping heterogeneity has not yet been systematically addressed.

Thus far, however, the relationship between the DEE molecular-level structure and the local geometry/organization and symmetry and size of the ingoing DEE components is not well understood. We suggest that a deeper understanding of how, e.g., the anion symmetry and size influence the HB network, which in turn affects the structural organization and macroscopic behavior, could lead to DEEs with enhanced physicochemical and electrochemical properties. We therefore here introduce molecular-level heterogeneity (MLH), defined mainly by the heterogeneous distribution of the HB network arising from variations in vital molecular properties such as size, charge, shape, and manifests as coordination and solvation distribution differences, and also in dynamic properties, such as ion diffusion rates. Hence the MLH causes non-uniform behaviour that ultimately impacts the overall performance of the DEEs.

To investigate how primarily anion symmetry and size affect the HB network and thus, in turn, cause MLH, we as a proof-of-concept study use three DEEs, with components defined above, at a 1[thin space (1/6-em)]:[thin space (1/6-em)]4 molar ratio by molecular dynamics (MD) simulations.

Methodology

xTB MD simulations

xTB simulations30,31 were implemented using CP2K32 with three-dimensional periodic boundary conditions using the geometry-, frequency-, noncovalent- and extended tight-binding (GFN-xTB) method.31 This semi-empirical quantum MD (xTB-MD) approach incorporates an augmented basis set with improved HB treatment and was specifically developed to yield accurate geometries, vibrational properties, and noncovalent interactions. It offers a favorable compromise between the efficiency of classical MD and the accuracy of DFT, enabling us to simulate the complex ion–solvent dynamics in DEEs.31

The initial configurations for the three DEEs were generated using CHAMPION,33 which relaxes the system using the conjugate gradient method to avoid atom overlaps and refines the geometry output. The 1[thin space (1/6-em)]:[thin space (1/6-em)]4 LiBF4/LiDFOB/LiBOB:NMA DEEs were composed of 486, 464, and 496 atoms, respectively (see SI for breakdown in molecules), and had cubic box side lengths of 16.7 Å, 17.3 Å and 16.7 Å, respectively, yielding initial densities within 0.98–1.1 g cm−3, consistent with previously reported values for similar systems. The total energy for the DEEs was obtained as:

 
E = Eel + Erep + Edisp + EXB(1)
where the total energy is represented by the electronic (Eel), atom-pairwise repulsion (Erep), dispersion (Edisp), and halogen-bonding (EXB) terms.31,34 For each DEE, an initial MD equilibration run (40 ps) was performed using an NPT ensemble with 1 fs time-steps at 363 K. The temperature was chosen to ensure that the DEE system remains liquid and to facilitate configurational sampling. Equilibration was assessed by monitoring the convergence of potential energy and the stability of temperature, volume, and RDFs, indicating that structural and thermodynamic properties had stabilized. The barostat set the pressure to 1 atm and a CSVR was used as a thermostat. The xTB method was used in both equilibration and production runs and once the energy and the density had converged, and evaluating the total energy as a function of simulation time, MD production runs (60 ps) were made in NVT ensemble with 1 fs time-steps and a Nosé–Hoover thermostat35 set to 363 K. Following this, the generated trajectories were analyzed using VMD,36 and the (local) structure studied by calculating the g(r) to generate radial distribution functions (RDFs) as:
 
image file: d5cp01828e-t1.tif(2)
where V represents the total volume, r is the distance between a pair of particles, p(r) is the average number of atom pairs found at a distance between r and r + dr and Npairs is the number of unique pairs of atoms where one atom is from each of two sets selections.37

The coordination number (CN) was subsequently calculated by integrating the g(r), over spherical coordinates, up to the first minimum of the RDF. Partial coordination numbers (pCNs) were similarly determined by integrating the partial RDFs, and also yield the cumulative distribution functions (CDFs) for specific pairs of atoms. The anion solvation numbers (SNs) were obtained by summing the pCN contributions of anion–solvent and anion–cation. In addition, we used Travis38 to generate spatial distribution functions (SDFs), which provide a three-dimensional representation.

To analyze the HB distribution, hydrogen bonds were first identified using geometric cutoffs of HBA–HBD distance ≤3.5 Å and HBA–H–D angle ≥150°. We then employed a Python script that from the MD simulation trajectories calculates the probability of finding an HB within a specific radius, based on HBD–HBA distances. The probability distribution function (PDF) was generated by binning these distances and normalizing the data to reflect the probability of HB formation. Additionally, we characterized the HB strength and spatial distribution using a contour plot through Travis,38 where the x-axis represents the HBD–HBA distances (e.g., O–H, F–H) and the y-axis the HB angles. To quantify “HB strength” and enable direct comparison across the three DEEs, we differentiated the cumulative distribution functions of these distance–angle bins to obtain a 2D probability density P(r,θ). All heat maps share a fixed color scale (vmin = 0; vmax = global maximum P over all DEEs), so identical color intensities correspond to identical absolute HB probabilities with warmer colors indicating more probable (stronger) HB geometries and cooler colors less frequent.

Classical MD simulations

MD simulations were performed to study dynamic properties and potential heterogeneity, as already roughly observed from the xTB simulations. The dynamic analysis is based on the calculation of self-diffusion coefficients from the mean-square displacements (MSDs).

The MD simulations were performed in LAMMPS,39 using the non-polarizable OPLSAA force field.40 To account for electronic screening effects and enhance the accuracy of inter-ionic interaction predictions, the partial charges of ions were scaled by a factor of 0.8.41,42 The simulation boxes contained 200 cations, 200 anions, and 800 NMA molecules. For each box, topology files, bonded and Lennard-Jones parameters, and initial configurations were generated using Packmol43 and fftool.44 The simulations consisted of an energy minimization step employing conjugate gradients, followed by equilibration in the microcanonical ensemble (NVE) for 500 ps, succeeded by another 1 ns in the isothermal–isobaric ensemble (NPT) at 1 atm, and an additional 1 ns in the canonical ensemble (NVT), all carried out at 363 K. These initial runs were designated as equilibration runs and were excluded from the subsequent analysis. Production runs were performed in the NVT ensemble for 14 ns at 363 K and 1 atm, using the average simulation box size from the NPT equilibration run. A Nosé–Hoover thermostat was used with a temperature damping of 100 fs and a pressure damping of 1000 fs. Electrostatic interactions were computed using the particle–particle–particle–mesh scheme (PPPM), and periodic boundary conditions were applied in all directions. MSDs were calculated using the compute msd subroutine in LAMMPS, which tracks the displacement of atoms over time and averages the squared distance between their current and reference positions.45 The self-diffusion coefficients were then calculated using a self-written script in MATLAB, determining the self-diffusion coefficient, D, via the Einstein relation:

 
image file: d5cp01828e-t2.tif(3)
where ri(t) is the position of the i-th particle at the time t, and d is the dimensionality of the system (e.g., d = 3 for three-dimensional diffusion). The brackets 〈·〉 denote the ensemble average, which is often approximated by a time average over the duration of the simulation. In practice, reasonable averages require sufficiently long simulation durations to ensure that the individual ion movements become uncorrelated and that the dynamics as monitored by the MSD reach the diffusive regime.46

Molecular-level heterogeneity

To quantify how primarily the HB network, but also other factors, causes MLH in the DEEs, an index is here proposed as a metric. This way, the MLH across the different DEEs can be compared, facilitating a deeper understanding. The index is designed to capture both structural and dynamic MLH, and is formulated as a weighted sum of selected features:
 
image file: d5cp01828e-t3.tif(4)
where fi represents individual molecular features, and wi are their corresponding weights, constrained such that image file: d5cp01828e-t4.tif. The features incorporated include the HB localization factor, the diffusion ratio (Danion/DNMA), the anion size, the anion SN, and symmetry elements. All features are normalized by their maximum value prior to PCA to prevent any single feature from dominating due to scale. This ensures fair weighting regardless of units. Features inversely proportional to heterogeneity (e.g., diffusion ratio, symmetry) are transformed viaimage file: d5cp01828e-t5.tif,aligning their interpretation with MHI directionality. Testing using z-score standardization showed preserved trends, with minor shifts in absolute weights. The final correlation-based weight refinement further corrects for any over- or under-weighting caused by scaling. The HB localization factor, derived from the Gini coefficient, which is a statistical measure of dispersion,47 is used here to quantify the HB distribution.

To compute the index, we employ a two-step hybrid approach combining Principal Component Analysis (PCA) and Correlation Analysis to determine the weights ωi. PCA is first used to assess which features contribute most to the variance across the DEEs. However, since PCA emphasizes statistical variance rather than physical relevance, we refine the weights through correlation analysis by computing the correlation between each feature and the initial MLH estimate. This hybrid method ensures that the final weights reflect both statistical significance and physicochemical meaning, aligning the index with heterogeneity-driving behaviors in the system.

The general trend of MHI values was preserved under alternative weighting schemes (e.g., uniform or PCA-only weights), but the hybrid approach improved numerical resolution and better reflected qualitative HB and mobility observations.

Results and discussion

We here present a detailed analysis of the local structure using the RDFs, the CNs, the anion SNs, and HB networks, derived from our xTB–MD trajectories to study various ion–ion and ion–solvent interactions (Fig. 2). Next, we analyze the dynamic behavior through MSDs and self-diffusion coefficients using classical MD simulations. Finally, we conclude with the MLH and connect it back to the structure and dynamics of the DEEs.
image file: d5cp01828e-f2.tif
Fig. 2 Schematic of the interactions: cation–anion, anion–solvent, and cation–solvent.

Local structure

We begin our analysis with the RDFs (Fig. 3), alongside the pCNs (Fig. 4) and SDFs (Fig. 5), to assess the local structure in the DEEs. We do this to scrutinize the roles of the different interactions present in the DEEs: cation–anion, cation–solvent, and anion–solvent interactions (Fig. 2).
image file: d5cp01828e-f3.tif
Fig. 3 RDFs for cation–anion interactions (a)–(c), cation–solvent interactions (d), (e), and anion–solvent interactions (f) for all the DEEs.

image file: d5cp01828e-f4.tif
Fig. 4 The three pCNs for each of the DEEs.

image file: d5cp01828e-f5.tif
Fig. 5 Snapshots after the production run of the: (a) LiBF4:NMA, (b) LiDFOB:NMA, and (c) LiBOB:NMA DEEs, respectively. Red is lithium cation, orange is anion, and blue is solvent.

From the RDFs (Fig. 3(a)–(c)), we observe that the lithium cation primarily interacts with the oxygen atoms of the DFOB and the BOB anions, with strong Li–O peaks at approximately 2 Å, indicating highly localized first cation coordination shells. The highest cation–anion pCN is observed in the LiBOB-based DEE (2.5), followed by the LiDFOB-based DEE and then the LiBF4 one (1.4) (Fig. 4). This trend follows the anion size and available coordination sites: BOB > DFOB > BF4.

The increased ion–ion interactions reflect the size of the anion.48

In more detail, in the LiDFOB-based DEE, the lithium cations tend to coordinate predominantly with the oxygen atoms (82%) rather than the fluorine atoms (18%), considering that a lack of preference would result in 66% and 33%. This selective interaction suggests that the presence of multiple coordination sites within DFOB facilitates more oxygen-centric solvation.

Turning to the cation–solvent interactions, the LiBF4-based DEE shows the highest lithium cation solvation by NMA (pCN = 2.7), with a sharp Li–O(NMA) peak at 1.9 Å (Fig. 3(e)), reinforcing the role of solvent structure. The cation–solvent solvation decreases in the order: LiBF4 > LiDFOB > LiBOB based DEEs, which hence inversely correlates with cation–anion coordination, and indeed the lower cation–solvent coordination in the LiBOB-based DEE we attribute to the stronger and dominant cation–anion interactions that limit the lithium cations’ ability to interact with NMA. Given that the total Li+ CN of ≈4.2 in the LiBF4 DEE, we note its agreement with the tetrahedral coordination found experimentally by Seo et al.49 in LiBF4·(AN)2 solvates comprising two acetonitrile nitrogen atoms and two BF4 fluorine atoms, indicating that Li+ favors a tetrahedral arrangement of its nearest neighbours.

Finally, for anion–solvent interactions, the B–N RDF (Fig. 3(f)) indicates that BF4, with its smaller size and higher symmetry, has the highest anion–solvent coordination (pCN = 3.0), followed by DFOB and then BOB (pCN = 0.7). The isotropic nature of BF4 facilitates uniform solvation, whereas the steric hindrance of BOB restricts solvent interactions, making it less solvated in comparison. The solvation trends overall suggest that the LiBOB-based DEE exhibits stronger ionic association, forming a distinct heterogeneous lithium–anion domain with the solvent at a longer distance, i.e., mostly outside the first cation solvation shell, leading to enhanced ion-pairing.

The main features in the B–B anion–anion RDF are somewhat closer and stronger, and also binodal for the LiBOB-based DEE at 4.1 Å and 5.0 Å (Fig. 6(a)).


image file: d5cp01828e-f6.tif
Fig. 6 RDFs of: (a) anion–anion, and (b) cation–cation correlations for the three DEEs.

This could suggest some kind of MH order, where these anions are localized together, something that is indeed visible in the snapshot (Fig. 5(c)).

Such binodal features are often associated with intermediate-range structuring or domain-like aggregation, consistent with MLH.

For the Li–Li cation–cation RDFs, there are no similar differences-all DEEs have the same main feature at approx. 3.5 Å (Fig. 6(b)). However, very noticeable is the much less ordered feature beyond this distance for the LiDFOB-based DEE, telling the story of the different cation solvation and cation–anion interactions in that particular DEE. These deviations in local ordering across salts support the manifestation of MLH, where local structural disorder varies based on anion symmetry and size.

And then, using the SDFs (Fig. 7) to visualize if there is a clear correlation between the choice of anion and the spatial distribution of ligands around them, we find that the LiBF4-based DEE, which exhibits the highest anion and cation solvation, also has the SDF showing the most “closed” and spatially uniform isosurface, thus following anion size and symmetry. Specifically, this enables more NMA nitrogen atoms to approach and interact with the anion, creating a rather homogenously structured DEE altogether.


image file: d5cp01828e-f7.tif
Fig. 7 SDFs for the: (a) LiBF4, (b) LiDFOB, and (c) LiBOB based DEEs showing lithium in purple and NMA in yellow around the anions.

In stark contrast, the SDFs for the DFOB and BOB-based DEEs display more localized and fragmented distributions, suggesting greater anisotropy and heterogeneity.

The cation coordination is more structured, especially for the LiBOB-based DEE, where distinct solvation “pockets” rather than a continuous shell are seen, due to the directional coordination preference introduced by the BOB anion.

HB analysis

In general, fluorine atoms tend to play a dominant role in HB formation in the DEEs, with F–H interactions being more prevalent than O–H interactions, due to their stronger electronegativity and spatial freedom. And indeed, in the LiDFOB-based DEE exhibits the most HBs due to its combination of fluorine and oxygen HBA sites, with pCN values of F–H = 1.3 and O–H = 0.8 (Fig. 8).
image file: d5cp01828e-f8.tif
Fig. 8 RDFs of the HBs in the three DEEs.

The distributions of the HB networks (Fig. 9) show that while both the LiDFOB-based and the LiBF4-based DEEs form similarly strong HBs, the latter shows slightly stronger HBs at shorter distances, corresponding to F–H distances of 1.5–2.0 Å (Fig. 9). Overall, the strongest HBs are found at shorter distances and larger angles (160–180°). In contrast, the BOB anion, being the largest anion with also directional preferences, sterically hinders HB formation and creates a localized HB network, as seen in the HB distribution, which shows a concentration at 6 Å (Fig. 9) and the HB strength, which shows the least intense heat map for its O–H HBs (Fig. 9). In line with these resuults, the HB lifetime analysis (Table S6) shows longer persistence for LiBOB and LiDFOB (O) compared to LiBF4 and LiDFOB (F), supporting the same trend.


image file: d5cp01828e-f9.tif
Fig. 9 Top: The HB strength heat map as a function of HBA–HBD distances r in pm and angles in degrees, with warmer colors indicating more probable (stronger) HB geometries. Bottom: The HB distributions of the three DEEs.

Dynamics

Turning swiftly to the dynamics of the DEEs, the total MSD curves (Fig. 10(a)) show an initial non-linear behavior, indicative of a ballistic regime followed by a diffusive motion.
image file: d5cp01828e-f10.tif
Fig. 10 (a) Total MSDs and (b) self-diffusion coefficients for the ions and the solvent.

The LiBF4-based DEE exhibits the steepest MSD slope (Fig. 10(a)), indicating the highest ion mobility, followed by LiDFOB and then LiBOB. This trend is also quantified in the inset of Fig. 11(a), where the mobility difference spans nearly a factor of two between LiBF4 and LiBOB electrolytes. This order is consistent with the calculated self-diffusion coefficients in Fig. 10(b), supporting that LiBF4 enables the fastest transport behavior.


image file: d5cp01828e-f11.tif
Fig. 11 (a) Radar plot comparing key molecular features: diffusion ratio, HB localization, anion size, and number of symmetry elements, and (b) bar chart of MHI for the LiBF4, LiDFOB, and LiBOB-based DEEs.

Fig. 10(b) further reveals distinct differences in species mobility not only across salts but also within each DEE. In the LiBF4-based DEE, the anion diffuses faster than the cation, with the solvent diffusing fastest overall. In contrast, the cation and anion in the LiDFOB-based and LiBOB-based DEE exhibit nearly identical self-diffusion coefficients, indicating stronger ion pairing and increasingly correlated motion. This narrowing of the cation–anion mobility gap from LiBF4 → LiDFOB → LiBOB DEEs reflects enhanced ionic association in the latter electrolytes. Similar lithium diffusivity has been reported for LiBF4 and LiTFSI in ADN-based electrolytes, closely matching the LiBF4-based DEE.50 Meanwhile, ethaline- and glycerol-based DEEs show slower lithium diffusion, comparable to or lower than those of LiDFOB and LiBOB in our study.51 In contrast, LiPF6-based electrolytes generally exhibit higher self-diffusion coefficients.52,53

Cross-comparing species within each electrolyte: in the LiBF4-based DEE, the trend observed is solvent > anion > cation; while in the LiDFOB-based and LiBOB-based DEEs, solvent > cation ≈ anion. The solvent mobility drops sharply between the LiBF4-based and the LiDFOB-based DEEs, after which it remains similarly suppressed in the LiBOB-based DEE. This shows that solvent motion becomes increasingly hindered as HB network localization and cation–anion ordering increase, suggesting a dynamically arrested environment in the most structured DEE.

The inverse correlation between localized HB networks (Fig. 8) and ion mobility (Fig. 10(b)), emphasizes the role of intermolecular interactions in determining not only extended structure but also transport properties.54,55

Molecular-level heterogeneity

Trying to get a comprehensive view of the MLH, we first gather some molecular features: HB localization factor, diffusion ratio, anion size, and number of symmetry elements (Fig. 11(a)), which do show quite some variability across the DEEs but do not allow for any clear discrimination between them.

Quantifying the structural and dynamic disorder within the DEEs by calculating an MHI (Fig. 11(b)), shows how the LiBOB- based DEE exhibits the largest MLH, due to its large anion size, strong cation–anion interactions, and localized HBs.

In contrast, the LiBF4-based DEE, with its higher ion solvation and more homogeneous HB distribution, shows the lowest MLH. Finally, the LiDFOB-based DEE, with a more moderate MHI suggests both spatial and dynamic heterogeneity, though to a lesser degree than the LiBOB-based DEE.

The physical meaning of the MHIs extend beyond relative comparison. Low MHIs (e.g., for LiBF4) reflect a more homogeneous molecular environment, associated with delocalized HB, relatively uniform diffusion behavior, and higher symmetry. In contrast, high MHIs (e.g., for LiBOB) indicate heterogeneous structuring and mobility, often arising from large, asymmetric anions and localized HB. Intermediate MHIs (e.g., for LiDFOB) capture partial heterogeneity in either structure or mobility, such as symmetry-breaking without full diffusion suppression. Thus, MHI serves as a semi-quantitative scale for molecular-level disorder, offering insights into dynamic uniformity or solvation–shell asymmetry, which is not apparent from single descriptors.

While the current MHI captures key aspects of MLH, it is worth noting that disorder can manifest itself in various forms. The index relies on five selected descriptors for interpretability, though future refinements could incorporate dynamic properties like HB lifetimes or local viscosity. PCA was limited to two components due to the small number of systems; expanding this would improve variance resolution. Lastly, although normalization by maxima reduces scaling bias, alternative standardization approaches may yield different weight distributions and deserve further investigation.

Conclusions

This research highlights how systematic variation in anion size and symmetry modulates solvation behavior, molecular interactions, transport dynamics, and molecular ordering in DEE-based electrolytes. The small, symmetric LiBF4-based DEE promotes isotropic coordination and uniform hydrogen bonding, enabling greater ion dissociation and solvent mobility. The LiBOB-based DEE, on the other hand, demonstrates the strongest cation–anion interactions, a localized HB network, and the highest MLH, with its ionic domains more segregated, possibly impacting both ion conduction pathways and perhaps even the electrochemical performance. The LiDFOB-based DEE falls, as expected, in between these two extremes, displaying moderate solvation, HB localization, and ionic mobility. Ultimately, these findings underscore the influence of anion structure on electrolyte properties, hydrogen-bond network arrangements, and MLH, providing valuable insights into the molecular-level interplay affecting electrolyte functionality. Such understanding can inform the rational design of future electrolytes tailored for improved ionic conduction and optimal electrochemical performance.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this study are provided in the SI. Supplementary information: This includes representative input files for both xTB (for production runs) and LAMMPS (for molecular dynamics equilibration and production runs). These files outline the key parameters and settings used to generate the simulation data discussed in the manuscript. All additional parameters, such as box sizes, system composition, etc, are described in the main text. Simulations were performed using the open-source software packages xTB (https://xtb-docs.readthedocs.io/) and LAMMPS (https://lammps.org/), and the generated trajectory data were used to compute structural and dynamic properties. The methods used for data analysis and post-processing are detailed in the main text, along with relevant equations and literature references, facilitating the simulation and analysis setup to be reproducible. Further instructions for running the input files and interpreting outputs. See DOI: https://doi.org/10.1039/d5cp01828e

Additional data or clarifications are available from the authors upon reasonable request.

Acknowledgements

The authors would like to acknowledge the financial support from the Swedish Research Council (grants #2020 – 03988 and #2021 – 00613) and VINNOVA/Batteries Sweden (BASE) (grant 2019-00064).

Notes and references

  1. J. Cabana, L. Monconduit, D. Larcher and M. R. Palacín, Beyond Intercalation-Based Li-Ion Batteries: The State of the Art and Challenges of Electrode Materials Reacting Through Conversion Reactions, Adv. Mater., 2010, 22, E170–E192 CrossRef CAS PubMed.
  2. A. Navarro- Suárez and P. Johansson, Perspective—Semi-Solid Electrolytes Based on Deep Eutectic Solvents: Opportunities and Future Directions, J. Electrochem. Soc., 2020, 167, 070511 CrossRef.
  3. A. Manthiram, X. Yu and S. Wang, Lithium battery chemistries enabled by solid-state electrolytes, Nat. Rev. Mater., 2017, 2, 1–16 Search PubMed.
  4. M. E. Di Pietro and A. Mele, Deep eutectics and analogues as electrolytes in batteries, J. Mol. Liq., 2021, 338, 116597 CrossRef CAS.
  5. A. Boisset, S. Menne, J. Jacquemin, A. Balducci and M. Anouti, Deep eutectic solvents based on N-methylacetamide and a lithium salt as suitable electrolytes for lithium-ion batteries, Phys. Chem. Chem. Phys., 2013, 15, 20054–20063 RSC.
  6. H. Liang, H. Li, Z. Wang, F. Wu, L. Chen and X. Huang, New binary room-temperature molten salt electrolyte based on urea and LiTFSI. The, J. Phys. Chem. B, 2001, 105, 9966–9969 CrossRef CAS.
  7. O. Geiculescu, D. DesMarteau, S. Creager, O. Haik, D. Hirshberg, Y. Shilina, E. Zinigrad, M. Levi, D. Aurbach and I. Halalay, Novel binary deep eutectic electrolytes for rechargeable Li-ion batteries based on mixtures of alkyl sulfonamides and lithium perfluoroalkylsulfonimide salts, J. Power Sources, 2016, 307, 519–525 CrossRef CAS.
  8. Z. Hu, F. Xian, Z. Guo, C. Lu, X. Du, X. Cheng, S. Zhang, S. Dong, G. Cui and L. Chen, Nonflammable nitrile deep eutectic electrolyte enables high-voltage lithium metal batteries, Chem. Mater., 2020, 32, 3405–3413 CrossRef CAS.
  9. L. J. B. M. Kollau, M. Vis, A. van den Bruinhorst, A. C. C. Esteves and R. Tuinier, Quantification of the liquid window of deep eutectic solvents, Chem. Commun., 2018, 54, 13351–13354 RSC.
  10. M. A. R. Martins, S. P. Pinho and J. A. P. Coutinho, Insights into the Nature of Eutectic and Deep Eutectic Mixtures, J. Solution Chem., 2019, 48, 962–982 CrossRef CAS.
  11. E. L. Smith, A. P. Abbott and K. S. Ryder, Deep Eutectic Solvents (DESs) and Their Applications, Chem. Rev., 2014, 114, 11060–11082 CrossRef CAS PubMed.
  12. J. Song, Y. Si, W. Guo, D. Wang and Y. Fu, Organosulfide-Based Deep Eutectic Electrolyte for Lithium Batteries, Angew. Chem., Int. Ed., 2021, 133, 9969–9973 CrossRef.
  13. Y. Hu, H. Li, X. Huang and L. Chen, Novel room temperature molten salt electrolyte based on LiTFSI and acetamide for lithium batteries, Electrochem. Commun., 2004, 6, 28–32 CrossRef CAS.
  14. S.-D. Han, J. L. Allen, E. Jónsson, P. Johansson, D. W. McOwen, P. D. Boyle and W. A. Henderson, Solvate structures and computational/spectroscopic characterization of lithium difluoro (oxalato) borate (LiDFOB) electrolytes, J. Phys. Chem. C, 2013, 117, 5521–5531 CrossRef CAS.
  15. C. W. Yang, Y. H. Ren, B. R. Wu and F. Wu, Formulation of a New Type of Electrolytes for LiNi1/3Co1/3Mn1/3O2 Cathodes Working in an Ultra-Low Temperature Range, Adv. Mater. Res., 2012, 455–456, 258–264 CAS.
  16. V. Aravindan, J. Gnanaraj, S. Madhavi and H. Liu, Lithium-Ion Conducting Electrolyte Salts for Lithium Batteries, Chem. – Eur. J., 2011, 17, 14326–14346 CrossRef CAS PubMed.
  17. J. Huang, L.-Z. Fan, B. Yu, T. Xing and W. Qiu, Density functional theory studies on the B-containing lithium salts, Ionics, 2010, 16, 509–513 CrossRef CAS.
  18. V. Aravindan and P. Vickraman, Effect of aging on the ionic conductivity of polyvinylidenefluoride–hexafluoropropylene (PVdF–HFP) membrane impregnated with different lithium salts, Indian J. Phys., 2012, 86, 341–344 CrossRef CAS.
  19. B. Bang and T. Yeu, Characterization of lithium tetrafluoroborate and N-methylacetamide complex as electric and ionic conductors, Korean J. Chem. Eng., 2016, 33, 1441–1446 CrossRef CAS.
  20. S. Kaur, A. Gupta and H. K. Kashyap, Nanoscale spatial heterogeneity in deep eutectic solvents, J. Phys. Chem. B, 2016, 120, 6712–6720 CrossRef CAS PubMed.
  21. Q. Li and H. Ardebili, Flexible thin-film battery based on solid-like ionic liquid-polymer electrolyte, J. Power Sources, 2016, 303, 17–21 CrossRef CAS.
  22. O. S. Hammond, D. T. Bowron and K. J. Edler, Liquid structure of the choline chlorideurea deep eutectic solvent (reline) from neutron diffraction and atomistic modelling, Green Chem., 2016, 18, 2736–2744 RSC.
  23. S. J. Bryant, A. J. Christofferson, T. L. Greaves, C. F. McConville, G. Bryant and A. Elbourne, Bulk and interfacial nanostructure and properties in deep eutectic solvents: Current perspectives and future directions, J. Colloid Interface Sci., 2022, 608, 2430–2454 CrossRef CAS PubMed.
  24. O. S. Hammond, D. T. Bowron and K. J. Edler, The effect of water upon deep eutectic solvent nanostructure: an unusual transition from ionic mixture to aqueous solution, Angew. Chem., Int. Ed., 2017, 129, 9914–9917 CrossRef.
  25. V. Alizadeh, D. Geller, F. Malberg, P. B. Sánchez, A. Padua and B. Kirchner, Strong microheterogeneity in novel deep eutectic solvents, ChemPhysChem, 2019, 20, 1786–1792 CrossRef CAS PubMed.
  26. D. O. Abranches and J. A. Coutinho, Everything you wanted to know about deep eutectic solvents but were afraid to be told, Annu. Rev. Chem. Biomol. Eng., 2023, 14, 141–163 CrossRef CAS PubMed.
  27. L. Percevault, A. Jani, T. Sohier, L. Noirez, L. Paquin, F. Gauffre and D. Morineau, Do deep eutectic solvents form uniform mixtures beyond molecular microheterogeneities?, J. Phys. Chem. B, 2020, 124, 9126–9135 CrossRef CAS PubMed.
  28. Y. Cui and D. G. Kuroda, Evidence of molecular heterogeneities in amide-based deep eutectic solvents, J. Phys. Chem. A, 2018, 122, 1185–1193 CrossRef CAS PubMed.
  29. S. Kaur, A. Malik and H. K. Kashyap, Anatomy of microscopic structure of ethaline deep eutectic solvent decoded through molecular dynamics simulations, J. Phys. Chem. B, 2019, 123, 8291–8299 CrossRef CAS PubMed.
  30. C. Bannwarth, S. Ehlert and S. Grimme, GFN2-xTB—An Accurate and Broadly Parametrized Self-Consistent Tight-Binding Quantum Chemical Method with Multipole Electrostatics and Density-Dependent Dispersion Contributions, J. Chem. Theory Comput., 2019, 15, 1652–1671 CrossRef CAS PubMed.
  31. S. Grimme, C. Bannwarth and P. Shushkov, A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z = 1–86), J. Chem. Theory Comput., 2017, 13, 1989–2009 CrossRef CAS PubMed.
  32. J. Hutter, M. Iannuzzi, F. Schiffmann and J. VandeVondele, cp2k atomistic simulations of condensed matter systems, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2014, 4, 15–25 CAS.
  33. R. Andersson, F. Årén, A. A. Franco and P. Johansson, CHAMPION: Chalmers hierarchical atomic, molecular, polymeric and ionic analysis toolkit, J. Comput. Chem., 2021, 42, 1632–1642 CrossRef CAS PubMed.
  34. S. Spicher and S. Grimme, Robust Atomistic Modeling of Materials, Organometallic, and Biochemical Systems, Angew. Chem., Int. Ed., 2020, 59, 15665–15673 CrossRef CAS PubMed.
  35. D. J. Evans and B. L. Holian, The Nose–Hoover thermostat, J. Chem. Phys., 1985, 83, 4069–4074 CrossRef CAS.
  36. W. Humphrey, A. Dalke and K. Schulten, VMD: Visual molecular dynamics, J. Mol. Graphics, 1996, 14, 33–38 CrossRef CAS PubMed.
  37. B. G. Levine, J. E. Stone and A. Kohlmeyer, Fast analysis of molecular dynamics trajectories with graphics processing units—Radial distribution function histogramming, J. Comput. Phys., 2011, 230, 3556–3569 CrossRef CAS PubMed.
  38. M. Brehm and B. Kirchner, TRAVIS-a free analyzer and visualizer for Monte Carlo and molecular dynamics trajectories, J. Chem. Inf. Model., 2011, 51(8), 2007–2023 CrossRef CAS PubMed.
  39. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19 CrossRef CAS.
  40. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, Development and Testing of the OPLS All-Atom Force Field on Conformational Energetics and Properties of Organic Liquids, J. Am. Chem. Soc., 1996, 118, 11225–11236 CrossRef CAS.
  41. E. S. C. Ferreira, I. V. Voroshylova, C. M. Pereira, M. Natália and D. S. Cordeiro, Improved Force Field Model for the Deep Eutectic Solvent Ethaline: Reliable Physicochemical Properties, J. Phys. Chem. B, 2016, 120(38), 10124–10137 CrossRef CAS PubMed.
  42. Z. Yu, P. E. Rudnicki, Z. Zhang, Z. Huang, H. Celik, S. T. Oyakhire, Y. Chen, X. Kong, S. C. Kim, X. Xiao, H. Wang, Y. Zheng, G. A. Kamat, M. S. Kim, S. F. Bent, J. Qin, Y. Cui and Z. Bao, Rational solvent molecule tuning for high-performance lithium metal battery electrolytes, Nat. Energy, 2022, 7, 94–106 CrossRef CAS.
  43. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, PACKMOL: A package for building initial configurations for molecular dynamics simulations, J. Comput. Chem., 2009, 30, 2157–2164 CrossRef PubMed.
  44. A. A. H. Pádua, Resolving dispersion and induction components for polarisable molecular simulations of ionic liquids, J. Chem. Phys., 2017, 146(20), 204501 CrossRef PubMed.
  45. S. Plimpton LAMMPS Documentation. https://www.lammps.org, 2024.
  46. D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, Elsevier, 2nd edn, 2002 Search PubMed.
  47. F. Cowell, Measuring Inequality, Oxford University Press, 2011 Search PubMed.
  48. M. Ue, A. Murakami and S. Nakamura, A Convenient Method to Estimate Ion Size for Electrolyte Materials Design, J. Electrochem. Soc., 2002, 149(10), A1385 CrossRef CAS.
  49. D. M. Seo, P. D. Boyle, J. L. Allen, S.-D. Han, E. Jónsson, P. Johansson and W. A. Henderson, Solvate Structures and Computational/Spectroscopic Characterization of LiBF4 Electrolytes, J. Phys. Chem. C, 2014, 118(32), 18377–18386 CrossRef CAS.
  50. O. Borodin and G. D. Smith, Mechanism of Ion Transport in Adiponitrile-Based Electrolytes: Molecular Dynamics Study of LiBF4 and LiTFSI in ADN, Phys. Chem. Chem. Phys., 2018, 20, 28294–28306 Search PubMed.
  51. Y. Ding, D. Sun and J. Xu, Lithium Ion Transport Properties in Deep Eutectic Solvents: A Molecular Dynamics Study of Ethaline and Glyceline-Based Electrolytes, Mol. Simul., 2021, 47, 1177–1186 Search PubMed.
  52. J. Landesfeind and H. A. Gasteiger, Temperature and Concentration Dependence of the Ionic Transport Properties of Lithium-Ion Battery Electrolytes, J. Electrochem. Soc., 2019, 166, A3079–A3097 CrossRef.
  53. A. B. Gunnarsdóttir, S. Vema, S. Menkin, L. E. Marbella and C. P. Grey, Investigating the Effect of a Fluoroethylene Carbonate Additive on Lithium Deposition and the Solid Electrolyte Interphase in Lithium Metal Batteries Using In Situ NMR Spectroscopy, J. Mater. Chem. A, 2020, 8, 14975–14992 RSC.
  54. U. L. Abbas, Q. Qiao, M. T. Nguyen, J. Shi and Q. Shao, Molecular dynamics simulations of heterogeneous hydrogen bond environment in hydrophobic deep eutectic solvents, AIChE J., 2022, 68, e17382 CrossRef CAS.
  55. M. T. Ong, O. Verners, E. W. Draeger, A. C. Van Duin, V. Lordi and J. E. Pask, Lithium ion solvation and diffusion in bulk organic electrolytes from first-principles and classical reactive molecular dynamics, J. Phys. Chem. B, 2015, 119, 1535–1545 CrossRef CAS PubMed.

This journal is © the Owner Societies 2025
Click here to see how this site uses Cookies. View our privacy policy here.