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
First published on 5th September 2025
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:
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.
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.
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.
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:
1 to 1
:
5. While there is no fixed ratio that guarantees the eutectic point, the optimal molar ratio for ionic conductivity was between 1
:
4 to 1
:
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:
4 molar ratio by molecular dynamics (MD) simulations.
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:
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) |
![]() | (2) |
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.
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:
![]() | (3) |
![]() | (4) |
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.
![]() | ||
Fig. 3 RDFs for cation–anion interactions (a)–(c), cation–solvent interactions (d), (e), and anion–solvent interactions (f) for all the DEEs. |
![]() | ||
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)).
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.
![]() | ||
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.
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.
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.
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
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.
Additional data or clarifications are available from the authors upon reasonable request.
This journal is © the Owner Societies 2025 |