Computational study of Li⁺ solvation structures in fluorinated ether, non-fluorinated ether, and organic carbonate-based electrolytes at low and high salt concentrations

Abstract

Understanding the solvation structure of electrolytes is crucial for optimizing the performance and stability of lithium-ion batteries. Novel electrolytes are essential for enhancing electrolyte structure and ensuring better integration with modern electrode systems. Herein, we report a novel weakly solvated ether electrolyte (WSEE) composed of a pure fluorinated ether solvent, which results in an anion-rich solvation structure even at a low salt concentration of 1 M. To explore this, we selected the advanced fluorinated solvent 2,2-difluoroethyl methyl ether (FEME) and compared it with dipropyl ether (DPE), ethylene carbonate (EC), and diethyl carbonate (DEC). The prepared electrolyte systems include DPE with 1 M, 1.8 M, and 4 M LiFSI; FEME with 1 M, 1.8 M, and 4 M LiFSI; and a 1 : 1 vol% EC/DEC mixture containing 1 M LiPF₆. In this work, we comprehensively investigate the Li⁺ solvation structures using molecular dynamics (MD) simulations and density functional theory (DFT) calculations. Our computational findings indicate the presence of large ion aggregates (AGGs) in each DPE- and FEME-based electrolyte, while SSIPs (68%) are the dominant species in the mixed EC/DEC electrolyte. Notably, the formation of large ion aggregates is more pronounced in FEME-based electrolytes. The dominant solvation structures in the ether-based electrolytes are the anion-rich complexes Li⁺(FSI⁻)₃(DPE)₁ and Li⁺(FSI⁻)₃(FEME)₁. We find that, similar to DPE, the FEME solvent also exhibits weak solvating power across all examined salt concentrations. More specifically, we find that FEME has weaker solvating power than DPE. This behavior is predicted by MD simulations, which indicate a strong preference for Li⁺ ions to coordinate with FSI⁻ anions within the primary solvation shell. We also observe that the number of unique solvation structures in the ether-based electrolytes increases with salt concentration, with FEME + LiFSI showing slightly more unique solvation structures than DPE + LiFSI. Furthermore, the quantum mechanical features of the Li⁺ solvation structures in DPE + 1.8 M LiFSI, FEME + 1.8 M LiFSI, and EC/DEC + 1 M LiPF₆ electrolytes are analyzed in detail using DFT calculations. We anticipate that this study will provide valuable insights into the Li⁺ solvation structures in DPE, FEME, and EC/DEC electrolytes, where the ether-based electrolytes exhibit closely similar properties.

---

1 Introduction

Rechargeable lithium-ion batteries (LIBs) are essential to modern energy storage, powering applications from consumer electronics to electric vehicles due to their high energy density, exceeding 300 Wh kg⁻¹,^1–7 long cycle life, and stability. Apart from electrode materials, choosing an appropriate electrolyte to facilitate lithium-ion transport between electrodes is a complex task that demands a deep understanding of the electrolyte structure.³ As the demand for longer-lasting and safer batteries grows, optimizing the electrolyte has become a critical research focus.⁸ To achieve this, frameworks are currently being developed to enable efficient searches for electrolyte materials.^9–11 The structure and dynamics of the electrolyte, particularly the solvation environment around Li⁺, play a vital role in determining key battery properties, including ion conductivity, electrochemical stability, and solid–electrolyte interphase (SEI) layer formation on the electrodes.^12–14 The transport mechanism of these Li⁺ ions within the electrolyte depends on their specific solvation structure, which is defined by the coordination of solvent molecules and anions around the Li⁺ ion. Therefore, a thorough understanding of the solvation structure is important for the development of improved electrolytes. Various analytical and computational techniques, including FTIR,^15–21 Raman,^22–24 NMR spectroscopy,^20–22,25 DFT calculations,^{19–21,26–28} and MD simulations^29–34 have been utilized to investigate solvation structures.

As mentioned above, the movement of Li⁺ ions through the electrolyte is recognized as a crucial factor influencing the rate at which energy is transferred to the electrodes.³⁵ According to the literature, Li⁺ ion transport occurs in two stages: first, the Li⁺ ions become surrounded by solvent molecules, and then these solvated ions migrate.³⁶ Gaining deeper insights into the solvation and transport behavior of Li⁺ ions can enable the development of improved electrolytes. Recent studies highlight the importance of the solvation structure of electrolytes and its impact on battery performance.^37,38 Specifically, weakly solvated ether electrolytes (WSEEs) have been shown to exhibit anion-rich solvation structures, which have attracted significant attention.^39,40 In contrast, carbonate-based electrolytes, widely used in commercial LIBs, often feature solvent-separated ion pairs (SSIPs) as the dominant solvation species due to their high solvation power.^36,41 This solvating power of a solvent is determined by the strength of the ion–dipole interaction between solvent molecules and Li⁺ ions. Chen et al. observed that the solvating power is governed by several factors, including the dipole moment and molecular orientation of the solvent, donor number, the extent of competition between solvents and anions in coordinating with Li⁺, and the dielectric constant.⁴² The study also showed that the competition between Li⁺–solvent and Li⁺–anion interactions largely determines the final Li⁺ solvation structures. However, Su et al. observed that the solvating power of solvents is primarily governed by their molecular structure, including steric hindrance and coordination ability, rather than by dielectric constant or donor number.⁴³ Earlier studies have further demonstrated that no single physical parameter, such as dielectric constant or dipole moment, can fully describe solvating power.⁴⁴ Chen et al. also noted that there is currently no clear consensus on how to define the solvating power of different electrolyte solvents.⁴²

Anion-rich solvation structures limit the interaction between free solvent molecules and Li⁺ ions, which helps suppress solvent decomposition and enhance electrolyte stability at both the anode and cathode interfaces.⁴⁵ These structures promote decomposition pathways dominated by the anion, leading to the formation of stable SEI layers enriched with inorganic components like LiF.⁴¹ As reported in the literature, these SEI layers play an important role in enhancing the performance of LIBs.^1,46 This enhancement is largely attributed to the SEI's ability to regulate Li⁺ ion migration at the electrode–electrolyte interface, which is governed by its composition and physicochemical properties.^47,48 Among its key features, the LiF-rich SEI layer functions as a robust protective shell on the electrode surface.⁴¹ Notably, an anion-rich solvation shell is known to facilitate the development of such LiF-rich SEI layers.⁴⁵ However, the strong coordination between Li⁺ and anions can raise the desolvation energy barrier, which may hinder Li⁺ transport—particularly in electrolytes with pronounced ion aggregation, such as those based on DPE or mixture of tetrahydrofuran (mixTHF).^41,45 Prior research has demonstrated that DPE/LiFSI-based electrolytes yield SEI layers with a high fluorine content (∼43%) and a significantly greater proportion of fluorinated species (∼22%) compared to other non-fluorinated ether systems like diethyl ether (DEE), 1,2-dimethoxyethane (DME), and diglyme (DIG).⁴⁵ Furthermore, DPE/LiFSI electrolytes are associated with enhanced cycling stability in lithium metal batteries.⁴⁵ These findings highlight the role of anion-driven interphase chemistry. Several studies have also focused on designing fluorinated electrolyte systems to facilitate the formation of weakly solvated structures and LiF-rich SEI layers.^49–54 Based on this insight, we have designed a novel fluorinated FEME ether-based electrolyte featuring an anion-abundant solvation structure to promote the formation of these desirable SEI layers.

In this study, we computationally investigate LiFSI-based non-aqueous electrolytes in fluorinated and non-fluorinated ether solvents, alongside highly soluble LiPF₆-based non-aqueous electrolytes in mixed carbonate solvents.⁵⁵ The investigation includes a series of solvents: DPE, FEME, EC, and DEC. For ether-based electrolytes, a broad range of salt concentrations (1 M, 1.8 M, and 4 M) is selected, while for the mixed carbonate electrolyte, only 1 M is studied. Our results demonstrate that DPE- and FEME-based electrolytes exhibit anion-rich solvation structures even at a low salt concentration of 1 M, with the results for the DPE electrolyte aligning with previous studies.⁴⁵ In contrast, the mixed carbonate electrolyte predominantly features homogeneously dispersed SSIPs.⁵⁶ We begin by analyzing the electronic properties of the salt and solvent molecules, including HOMO/LUMO distributions and electrostatic potential (ESP) maps, to evaluate their chemical stability within the electrolyte. Next, fluorinated, non-fluorinated, and mixed carbonate-based electrolyte systems are modeled using molecular dynamics (MD) simulations to derive key properties of interest. Radial distribution functions (RDFs) and coordination numbers (CNs) are then calculated to provide insight into the solvation structure of these systems. The solvation structures of lithium ions are examined in detail, with a focus on the composition of the primary solvation shell, including the number of solvent molecules and anions. Additionally, the presence of solvent-separated ion pairs (SSIPs), contact ion pairs (CIPs), and aggregated species (AGGs) is analyzed for each electrolyte. Electronic characteristics of Li⁺ solvation structures are further explored using charge density difference analysis, Bader charge calculations, electrostatic potential maps, binding energies, and HOMO/LUMO distributions derived from density functional theory (DFT) calculations. This comprehensive approach provides valuable insights into the solvation behavior of fluorinated, non-fluorinated ether, and mixed carbonate-based electrolytes, contributing to the design of advanced electrolytes for lithium-ion battery applications.

2 Computational methods and details

2.1 MD simulations

In this study, non-fluorinated ether, fluorinated ether, and mixed organic carbonate solvents were considered. The selected solvents were dipropyl ether (DPE), 2,2-difluoroethyl methyl ether (FEME), ethylene carbonate (EC), and diethyl carbonate (DEC). Herein, seven types of electrolytes were studied: DPE + 1 M LiFSI, DPE + 1.8 M LiFSI, DPE + 4 M LiFSI, FEME + 1 M LiFSI, FEME + 1.8 M LiFSI, FEME + 4 M LiFSI, and 1 : 1 vol% EC/DEC + 1 M LiPF₆. The Li⁺ solvation structures in these electrolytes were investigated using classical molecular dynamics (MD) and density functional theory (DFT) simulations.

In each simulation box, lithium cations Li⁺ and FSI⁻ anions were randomly distributed among 542 solvent molecules for the DPE + 1 M LiFSI, DPE + 1.8 M LiFSI, and DPE + 4 M LiFSI systems. For the FEME + 1 M LiFSI, FEME + 1.8 M LiFSI, and FEME + 4 M LiFSI systems, the ions were randomly placed among 787 solvent molecules. Similarly, the 1 : 1 vol% EC/DEC + 1 M LiPF₆ system contained 875 solvent molecules with Li⁺ and PF₆⁻ ions arranged in a non-uniform manner. Table 1 shows the density and molar mass of the solvent molecules.^66,67 In Table 2, the number of salt and solvent molecules in the electrolytes for any specific concentration was calculated (Note 1, ESI†).

Table 1 Density and molar mass of solvent molecules

Solvent	Density (g m^−3)	Molar mass (g mol^−1)	#Atoms
DPE	736 × 10^3	102.177	21
FEME	1004 × 10^3	96.076	12
EC	132 × 10^4	88.06	10
DEC	975 × 10^3	118.13	18

---

Table 2 Number of salt and solvent molecules in each electrolyte

Electrolyte	Concentration (M)	#DPE/#FEME/#EC	#DEC	#LiFSI/#LiPF6	#Atoms
DPE	1	542	—	75	12 132
DPE	1.8	542	—	135	12 732
DPE	4	542	—	301	14 392
FEME	1	787	—	75	10 194
FEME	1.8	787	—	135	10 794
FEME	4	787	—	301	12 454
1 : 1 EC/DEC	1	564	311	75	11 838

---

The initial configurations of these seven electrolyte systems were modeled using the PACKMOL⁶⁸ package by randomly placing the solvent molecules, FSI⁻, PF₆⁻, and Li⁺ in a 5 × 5 × 5 nm³ cubic simulation box. Atomic and ionic interactions in the DPE-LiFSI, FEME-LiFSI, and EC-DEC-LiPF₆ systems were described by the OPLS-AA (optimized potentials for liquid simulations all atom)^69,70 and ionic liquids force field. The bonded and non-bonded parameters of the OPLS-AA force field for the DPE, FEME, EC, and DEC solvent molecule atoms were obtained from LigParGen,^71–73 while the force field parameters for FSI⁻, PF₆⁻, and Li⁺ ions were obtained from a database of several ionic liquids.^74,75 This OPLS-AA force field has been extensively validated for modeling lithium-ion battery electrolytes and offers a favorable balance between accuracy and computational efficiency.

The functional form of the OPLS force field is defined by a set of potential functions³⁶ in eqn (1)–(7), which include:

E_total(r^N) = E_bonded + E_nonbonded (1)

E_bonded = E_bonds + E_angles + E_dihedrals + E_impropers (2)

where,

In MD simulations, interactions within molecules are divided into bonded and non-bonded interactions. Bonded interactions, which include bonds, angles, dihedrals, and impropers, are modeled using harmonic functions. Non-bonded interactions include van der Waals forces and Coulombic forces, which describe the behavior between atoms that are not directly bonded. The dihedral term captures the torsional motion of four consecutively bonded atoms, and the improper term describes the torsional motion of three atoms arranged around a central fourth atom.

All MD simulations were performed using the LAMMPS⁷⁶https://lammps.org open-source software (version 23 Jun 2022). Lennard-Jones and Coulombic force interactions were cut off at a distance of 1.2 nm. Coulombic forces beyond the cutoff were computed using the particle–particle particle–mesh (PPPM)⁷⁷ method to account for long-range electrostatic interactions with a relative error in forces of 1 × 10⁻⁵. Periodicity was applied in all the x, y, and z dimensions of the cubic simulation box. The equilibration procedure and production run are outlined as follows.⁴⁵ First, the prepared systems from the PACKMOL software were minimized using the steepest descent (SD) method with a convergence criterion of 1000 kcal mol⁻¹ Å⁻¹, followed by conjugate gradient (CG) minimization with a convergence criterion of 10 kcal mol⁻¹ Å⁻¹. Minimization algorithms were used to reduce the system's energy and prevent particle overlap. The systems were then equilibrated at a temperature of 298.15 K and a pressure of 1 atm in the isobaric–isothermal (NPT) ensemble using a time step of 1 fs for 2 ns to stabilize the potential energy and density of the systems. During equilibration, bond constraints were applied to specified bond lengths in the simulation using the SHAKE algorithm.⁷⁸ The temperature and pressure were regulated by the Nosé–Hoover thermostat and barostat,^79–81 with time constants set to produce characteristic fluctuations over 100 and 1000 time steps, respectively. Next, the equilibrated systems were heated to 500.15 K for 2 ns and then gradually cooled to 298.15 K over four steps, spanning 3 ns. Finally, production runs were performed in the canonical (NVT) ensemble at 298.15 K for 5 ns using a time step of 1 fs, from which the properties of interest were derived. The Nosé–Hoover thermostat was used in the NVT ensemble. All MD simulations and DFT calculations were carried out using our HPC cluster Wulver at NJIT and the Expanse supercomputing cluster at SDSC.

2.2 DFT calculations

In this research, all first-principles calculations in DFT were performed using the Vienna ab initio simulation package (VASP) software.^82,83 The core–valence electron interactions were treated using the projector augmented wave method or PAW pseudopotentials.^84,85 The commonly used Perdew–Burke–Ernzerhof (PBE) functional⁸⁶ under the Generalized Gradient Approximation (GGA) was used to model exchange and correlation interactions between electrons.

To perform geometry optimization in VASP, the atomic coordinates were allowed to change while keeping the shape and volume of the cell constant. Gaussian smearing was employed, with the smearing value set to 0.05 eV. The self-consistent field energy convergence was set to 1 × 10⁻⁶ eV, and the ionic force convergence tolerance was set to 0.02 eV Å⁻¹. An energy cutoff for the plane-wave basis set was specified at 520 eV. The Brillouin zone of the supercell was sampled using Γ-centered k-point grids (KPOINTS) of 3 × 3 × 3 with a k-mesh density of 0.03. Since the PBE functional provides a poor description of dispersion forces, the zero-damping DFT-D3 method of Grimme⁸⁷ was implemented to more accurately calculate the energy of the system. For calculations involving individual molecules, the size of the periodic models was set to 10 × 10 × 10 ų, while for clusters, the size was set to 15 × 15 × 15 ų and 20 × 20 × 20 ų. This setup ensures sufficient vacuum distance without significantly increasing the computational cost. In this study, all DFT calculations were conducted in vacuum. Both the geometry optimization and single point energy calculation, also called the self-consistent field (SCF) calculation were performed using non-spin-polarized calculations. The electronic structure information was obtained from SCF calculations performed on the optimized structures. This includes the ESP map, HOMO/LUMO distribution, Bader charge analysis, and charge density difference (CDD) of all the studied molecules and solvation structures.^{45,59,88–93} The binding energies were also calculated using DFT in vacuum with cluster models.^45,94,95 For structure visualization, VESTA⁹⁶ and Maestro were utilized, while the VASPKIT⁹⁷ package was employed for post-processing the wave functions (WAVECAR) to generate the HOMO and LUMO distributions. The ESP maps, HOMO/LUMO diagrams, and CDD were also plotted using VESTA software.

3 Results and discussion

3.1 Salts, solvents, and electrolyte systems

A variety of solvents, including non-fluorinated ether (DPE), fluorinated ether (FEME), and carbonate solvents (EC and DEC) were selected to investigate the Li⁺ solvation structures and coordination in the electrolyte (Fig. 1 and Table 1). All MD production runs were conducted at 25 °C. The ether-based electrolytes used LiFSI, while the mixed carbonate-based electrolyte was prepared with LiPF₆ salt. To systematically compare solvation structures across different solvent environments, we selected LiFSI for both DPE and FEME ether-based electrolytes due to its tendency to promote anion-rich solvation structures even at relatively low concentrations. The combination of fluorinated FEME with LiFSI contributes to the formation of a robust LiF-rich SEI layer, which is critical for interfacial stability.⁴¹ Prior studies have also reported anion-rich solvation in similar WSEE systems, such as DPE + 1.8 M LiFSI, DEE + 1.8 M LiFSI, and FDMB + 1 M LiFSI.^45,98,99 In contrast, LiPF₆ was selected for the carbonate-based EC/DEC system due to its widespread use in commercial LIBs and its formation of SSIPs.⁵⁶ This makes it a meaningful benchmark for comparison with our ether-based systems. Using LiFSI salt allowed us to compare the anion-rich solvation structure of the FEME electrolyte with the SSIP-dominated carbonate-based electrolyte (LiPF₆ in EC/DEC), while also enabling a controlled comparison with non-fluorinated DPE to highlight the enhanced ion aggregation promoted by the FEME + LiFSI system.

Fig. 1 2D and 3D structures with ESP maps of the studied salts and solvents. The red and blue regions represent areas of high electron density (negative charge) and low electron density (positive charge), respectively. The magenta, green, orange, red, gray/brown, white/light pink, yellow, and blue spheres represent Li, F, P, O, C, H, S, and N atoms, respectively. For all ESP maps, the isosurface level was set at 10 eV. The isosurfaces of ESP maps were visualized using VESTA.

DFT calculations in VASP were used to simulate the electrostatic potential (ESP) maps, highest occupied molecular orbital (HOMO), and lowest unoccupied molecular orbital (LUMO) of the salt/solvent molecules. Fig. 1 shows regions of high electron density (negative charge) and low electron density (positive charge) in the ESP maps, with values ranging from 0 to 1.24 e bohr⁻³ (DPE, EC, DEC) and from 0 to 2.08 e bohr⁻³ (FEME, LiFSI, LiPF₆). The isosurfaces and energy levels of the HOMO/LUMO distributions for DPE, FEME, EC, DEC, LiFSI, and LiPF₆ molecules are depicted in Fig. 2. The higher the HOMO, the easier it is for the molecule to donate electrons, while a lower LUMO indicates it can accept electrons more easily. The energy difference between the HOMO and LUMO is referred to as the energy band gap (LUMO–HOMO). A narrower band gap usually corresponds to greater chemical reactivity and lower stability, while a wider band gap suggests less reactivity and greater stability. Herein, the lower HOMO (−6.2697 eV) and LUMO (−0.9929 eV) of LiFSI compared to DPE and FEME indicate it will decompose first during the charge/discharge cycle.⁵⁹ Likewise, the lower HOMO (−8.0773 eV) and LUMO (−1.1581 eV) of LiPF₆ relative to EC and DEC suggest LiPF₆ will decompose before the solvents and earlier than LiFSI. According to the literature, slight variations in the HOMO/LUMO energy values of the same molecule can result from using different calculation methods (PBE, GGA, B3LYP functional) and software packages such as VASP and Gaussian.^{57–65,100–107} The DFT validations are shown in Fig. 3. Our HOMO/LUMO energy values of the salts and solvents are quite similar to those obtained using the PBE functional and VASP.

Fig. 2 Comparison of HOMO/LUMO energy levels of the studied salts and solvents. The positive and negative phase of HOMO and LUMO are depicted in yellow and cyan colors, respectively. Yellow and cyan indicate the different signs of the isosurface of the wave function, and their sizes indicate its amplitude. For all HOMO/LUMO diagrams, the isosurface level was set between 8 × 10⁻⁹ and 8 × 10⁻⁸ e bohr⁻³. The isosurfaces of HOMO and LUMO were visualized using VESTA.

Fig. 3 DFT validation of HOMO/LUMO energy levels of the studied salts and solvents. This work, Wu et al.,⁵⁷ Jiang et al. (2023a),⁵⁸ and He et al.⁵⁹ used the PBE functional in VASP, while Vangapally et al.,⁶⁰ Hou et al.,⁶¹ Jiang et al. (2023b),⁶² Ma et al.,⁶³ Wang et al.,⁶⁴ and Chen et al.⁶⁵ used the B3LYP functional in Gaussian software.

The electrolyte systems were modeled using MD simulations in LAMMPS. Each system followed a process of minimization, equilibration, melting, quenching, and a 5 ns production run, from which key properties such as radial distribution functions (RDF), coordination numbers (CN), and solvation structures were obtained. Fig. 5 shows the 3D simulation box for the 1.8 M and 1 M DPE-LiFSI, FEME-LiFSI, and EC-DEC-LiPF₆ electrolyte systems after the production run. The ESI† includes OVITO-generated GIF showing the 5 ns production run trajectories for these electrolyte systems.¹⁰⁸ The 3D simulation boxes of all the remaining systems, including DPE-LiFSI and FEME-LiFSI at 1 M and 4 M concentrations are depicted in Fig. S1 (ESI†). Additionally, Fig. S2–S5 (ESI†) provide the volume, density, pressure, and temperature of all the equilibrated systems. During equilibration, the pressure and temperature fluctuated around 1 atm and 298.15 K, respectively. Larger systems with more atoms (4 M) exhibited fewer fluctuations compared to smaller systems (1 M and 1.8 M). The initial approximate concentration (before NPT equilibration), the actual concentration after NPT equilibration, and the system density during the NPT equilibration process are shown in Fig. 4. The density of our equilibrated system, 1 M LiPF₆ in a 1 : 1 EC/DEC mixture is 1.240 g cm⁻³ at 25 °C. This value is in good agreement with experimental data reported in the literature: Anton Paar reports 1.242 g cm⁻³ at 20 °C, and Sigma-Aldrich lists 1.26 g cm⁻³ at 25 °C. Additionally, published values include 1.225 g cm⁻³ at 25 °C from Lundgren et al., 1.23 g cm⁻³ from Lee et al., and 1.25 g cm⁻³ from Dougassa et al., all at 25 °C.^109–112

Fig. 4 (a) Initial approximate concentration before NPT equilibration and post-NPT concentration of the equilibrated system. (b) Density of the system during NPT equilibration.

Fig. 5 Snapshots of the simulation box obtained from MD simulations at 5 ns of the production run: (a) DPE + 1.8 M LiFSI, (b) FEME + 1.8 M LiFSI, and (c) 1 : 1 EC/DEC + 1 M LiPF₆.

3.2 Radial distribution function and coordination number

The radial distribution function (RDF) and coordination number (CN) of the electrolytes were calculated during MD simulations to study the electrolyte structure and the coordination environment of Li⁺. The RDF and CN for Li⁺–O(DPE), Li⁺–O(FEME), Li⁺–O(EC), Li⁺–O(DEC), Li⁺–O(FSI⁻), Li⁺–N(FSI⁻), Li⁺–P(PF₆⁻), Li⁺–F(PF₆⁻), and Li⁺–Li⁺ were calculated by averaging across 1 million configurations of each electrolyte system during the last 1 ns of the 5 ns production run. These results are shown in Fig. 6 and 7 (Fig. S6 and S7, ESI†), with the corresponding numerical data provided in Table 3. For more detailed information, Fig. S8–S14 (ESI†) also provide the time-averaged RDF and CN plots for each 1 ns interval of the 5 ns production run. Since the RDF and CN remained consistent across each 1 ns interval, a total production run of 5 ns was selected.

Fig. 6 (a) RDFs of Li⁺–O(DPE/FEME/EC/DEC) and corresponding (b) coordination numbers as a function of distance. (c) RDFs of Li⁺–O(FSI⁻), Li⁺–F(PF₆⁻), and corresponding (d) coordination numbers as a function of distance.

Fig. 7 (a) RDFs of Li⁺–Li⁺ and corresponding (b) coordination numbers as a function of distance. (c) RDFs of Li⁺–N(FSI⁻), Li⁺–P(PF₆⁻), and corresponding (d) coordination numbers as a function of distance.

Table 3 RDF and coordination in ether and carbonate-based electrolytes. Cutoff for RDF is the maximum peak; for CN, it is the first minimum after the first peak in RDF (first solvation shell)

Electrolyte	M	Pair	RDF	Cutoff (Å)	CN	Cutoff (Å)
DPE	1 M	Li–O(DPE)	17.25	1.97	0.77	2.99
		Li–O(FSI)	106.15	2.09	3.82	3.00
		Li–Li	9.84	5.83	3.95	6.93
		Li–N(FSI)	23.25	4.46	3.47	5.14
DPE	1.8 M	Li–O(DPE)	17.44	1.97	0.80	3.00
		Li–O(FSI)	61.99	2.09	3.81	3.00
		Li–Li	5.50	5.43	3.86	7.01
		Li–N(FSI)	13.48	4.50	3.43	5.26
DPE	4 M	Li–O(DPE)	20.64	1.97	0.77	3.00
		Li–O(FSI)	33.45	2.08	3.75	3.00
		Li–Li	2.56	6.00	3.94	7.05
		Li–N(FSI)	7.30	4.45	3.33	5.23
FEME	1 M	Li–O(FEME)	7.70	2.02	0.72	2.98
		Li–O(FSI)	103.43	2.09	3.92	2.97
		Li–Li	9.57	5.23	5.34	7.27
		Li–N(FSI)	28.48	4.49	3.90	5.56
FEME	1.8 M	Li–O(FEME)	7.98	2.02	0.74	3.00
		Li–O(FSI)	62.14	2.09	3.96	3.00
		Li–Li	5.76	5.94	5.03	7.11
		Li–N(FSI)	17.15	4.49	3.86	5.31
FEME	4 M	Li–O(FEME)	9.69	2.02	0.79	3.00
		Li–O(FSI)	31.84	2.09	3.91	3.00
		Li–Li	2.85	5.93	4.96	7.22
		Li–N(FSI)	7.93	4.49	3.63	5.19
1 : 1 EC/DEC	1 M	Li–O(EC)	8.41	2.11	0.90	3.49
		Li–O(DEC)	140.33	1.99	3.58	2.99
		Li–F(PF6)	6.86	2.06	0.56	2.91
		Li–Li	0.26	5.93	0.09	7.32
		Li–P(PF6)	6.54	3.57	0.38	4.79

---

In our simulations, weak solvating power was identified in ether solvents (DPE and FEME), as indicated by their low RDF peaks, low coordination numbers, and weak binding energies with Li⁺. Conversely, in the mixed carbonate systems, EC and DEC exhibited stronger coordination with Li⁺, evident from sharper RDF peaks, higher CNs, and greater binding energies, suggesting stronger solvating power. This contrast highlights the impact of solvent type—ether versus carbonate—on solvation behavior and ion coordination. Such findings are crucial for tailoring electrolytes to achieve desired interfacial properties in LIBs. These results show that solvating power is not determined by a single parameter but arises from a combination of molecular properties and competitive ion interactions within the solvation shell.^42–44 In this work, two carbonate solvents in a mixed-solvent electrolyte, with a solvating power order of EC < DEC (computationally), and two ether solvents in a single-solvent electrolyte, with a solvating power order of FEME < DPE were investigated.

The RDF and CN were calculated using eqn (8) and (9), where n(r) is the average number of particles in the spherical shell, r is the interatomic separation distance, and ρ is the particle density in the system. As reported in the literature, the cutoff or threshold value of the pairwise distance (r) for the RDF corresponds to the distance at the maximum peak, while for the CN, it is the first minimum following the first peak in the RDF, representing the first solvation shell.⁵⁶

In this study, three different salt concentrations (1 M, 1.8 M, and 4 M) were used for both DPE and FEME electrolytes, whereas only a 1 M concentration was applied to the EC/DEC electrolyte. The Li⁺–solvent and Li⁺–anion pairwise interactions follow similar trends in each ether electrolyte but exhibit opposite trends in the carbonate electrolyte (Fig. 6, 7 and Table 3). In the mixed EC/DEC electrolyte, DEC solvents exhibit stronger coordination with Li⁺ compared to EC, as indicated by the sharp RDF peaks of Li⁺–O(DEC) (Fig. 6a and b).^20,50,51,56 Conversely, the minimal Li⁺–F(PF₆⁻) RDF peaks suggest that the PF₆⁻ anion has little influence on Li⁺ solvation structures (Fig. 6c and d). On the other hand, Fig. 6 shows that in both DPE and FEME electrolytes, due to their weak solvating power as indicated by the negligible RDF peaks of Li⁺–O(DPE) and Li⁺–O(FEME), DPE and FEME solvents do not strongly interact with the cations. This allows Li⁺ to form strong interactions with the FSI⁻ anion in the primary solvation shell (r ≈ 2.09 Å), as shown by the pronounced RDF peaks of Li⁺–O(FSI⁻).^{45,102,113,114} From Fig. 6a, c, and Table 3, in both DPE and FEME electrolytes, as salt concentration increases, the maximum RDF peaks for Li⁺–O(DPE) and Li⁺–O(FEME) increase, while those for Li⁺–O(FSI⁻) decrease. It is found that the maximum RDF peaks of Li⁺–O(FSI⁻) are consistently higher than those of Li⁺–O(DPE) and Li⁺–O(FEME), which further confirms the weak solvating power of DPE and FEME solvents across varying salt concentrations. Additionally, the lower RDF peaks and coordination number of Li⁺–O(FEME) compared to Li⁺–O(DPE) suggest that FEME has weaker solvating power than DPE. This observation is also supported by the coordination number of Li⁺–O(FSI⁻), which is slightly higher in FEME (ranging from 3.91 to 3.96) than in DPE (ranging from 3.75 to 3.82) electrolytes (Table 3). The cutoffs for the sharp RDF peaks of Li⁺–O(DPE), Li⁺–O(FEME), and Li⁺–O(FSI⁻) are observed at r = 1.97 Å, r = 2.02 Å, and r = 2.09 Å, respectively, and remain unchanged across different salt concentrations.

From Fig. 6c, 7a and c, the RDFs of Li⁺–O(FSI⁻), Li⁺–F(PF₆⁻), Li⁺–N(FSI⁻), Li⁺–P(PF₆⁻), and Li⁺–Li⁺ suggest the formation of ion clusters in all the studied electrolytes. Fig. 7a shows that shorter Li⁺–Li⁺ interaction distances are observed at r = 5.83 Å, 5.43 Å, 6.00 Å (1st peak) and 7.85 Å, 8.15 Å, 8.17 Å (2nd peak) for the DPE electrolytes, and r = 5.23 Å, 5.94 Å, 5.93 Å (1st peak) and 8.04 Å, 8.08 Å, 8.09 Å (2nd peak) for the FEME electrolytes, corresponding to 1 M, 1.8 M, and 4 M concentrations, respectively. Additionally, in these ether electrolytes, shorter Li⁺–N(FSI⁻) interaction distances are also observed at r ≈ 2.20 Å and r ≈ 4 Å, where two consecutive RDF peaks are identified around 4 Å (Fig. 7c). These shorter interaction distances of Li⁺–Li⁺ and Li⁺–N(FSI⁻) indicate the presence of large ion aggregates (AGGs) composed of multiple Li⁺ and FSI⁻ ions.^45,115 Across all the salt concentrations (1 M, 1.8 M, and 4 M) in DPE and FEME ether electrolytes, the Li⁺ and FSI⁻ ions form large ion aggregates through bridging coordination, where the FSI⁻ anions coordinate with multiple Li⁺ cations via their O atoms (Fig. 5 and Fig. S1, ESI†). This aggregation behavior is further supported by the long-range Li⁺–O(FSI⁻) interactions in the secondary solvation shell, observed at r ≈ 4.38 Å in DPE and 4.36 Å in FEME electrolytes (Fig. 6c). In contrast, the EC/DEC carbonate electrolyte exhibits a homogeneous distribution of Li⁺, PF₆⁻, and EC/DEC solvent molecules (Fig. 5). The Li⁺ cations in the carbonate electrolyte are more widely separated as they are strongly coordinated by solvent molecules, which leads to the formation of solvent separated ion pairs (SSIPs).^{26,41,56,116–119} This is confirmed by the lower RDF peaks and coordination numbers of Li⁺–Li⁺ interactions in the primary solvation shell (Fig. 7a and b). Similarly, the lower RDF peaks and CN of Li⁺–F(PF₆⁻) and Li⁺–P(PF₆⁻) interactions also confirm the significant presence of SSIPs in this carbonate electrolyte (Fig. 6c, d, 7c, and d). In this study, for all seven electrolytes, the cutoff value used to calculate the coordination number of Li⁺–O(solvent) and Li⁺–O(anion) is approximately 3 Å (3.49 Å for Li⁺–O(EC)), and for Li⁺–Li⁺, it is approximately 7 Å (Table 3). Within this cutoff, the CN of Li⁺–O(solvent) is lower in FEME (0.72, 0.74, 0.79) compared to DPE (0.77, 0.80, 0.77), while the CN of Li⁺–O(FSI⁻) is higher in FEME (3.92, 3.96, 3.91) than in DPE (3.82, 3.81, 3.75) at 1 M, 1.8 M, and 4 M concentrations, respectively. Similarly, the CN of Li⁺–Li⁺ interactions is higher in FEME (5.34, 5.03, 4.96) than in DPE (3.95, 3.86, 3.94), whereas it is significantly lower in the EC/DEC system (0.09). In all electrolytes, the RDF and CN of Li⁺–Li⁺ interactions remain zero up to r = 3 Å (Fig. 7a and b). The atomic coordination numbers for all electrolytes are also compared in Fig. 8. Our findings reveal that DPE and FEME electrolytes are primarily composed of AGGs, while SSIPs are more prevalent in the EC/DEC electrolyte. Inside the first solvation shell, the higher RDF peaks and CN values of Li⁺–N(FSI⁻) in FEME electrolytes also highlight that Li⁺–FSI⁻ pairs exhibit stronger aggregation in FEME electrolytes than in DPE electrolytes across all the salt concentrations (Fig. 5, 7c, d, Fig. S1 and Table 3, ESI†). This is further supported by the higher coordination numbers of Li⁺–Li⁺ and Li⁺–O(FSI⁻) interactions in the primary solvation shell of FEME electrolytes compared to all other electrolytes (Table 3). Our RDF results and solvation structures for DPE + 1.8 M LiFSI using the OPLS-AA^69,70 force field are quite similar to those reported by Li et al. using the OPLS-2005 force field from Schrödinger.^45,120 The properties of the 1 : 1 EC/DEC + 1 M LiPF₆ electrolyte also closely match those reported by Gullbrekken et al., with both studies using the OPLS-AA force field.⁵⁶

Fig. 8 Atomic coordination numbers of the electrolytes. The CN of Li⁺–Li⁺ and Li⁺–F(PF₆⁻) in mixed EC/DEC carbonate electrolyte is assumed to be the same for both EC and DEC solvents.

3.3 Solvation structure

In this work, all lithium-ion solvation structures were obtained from MD snapshots. These solvation structures significantly influence ion transport, stability, and chemical reactivity in electrolytes, which are critical for the design of lithium-ion batteries. The solvation structures in DPE- and FEME-based electrolytes are categorized into three groups: AGG-1 (containing one or two FSI⁻ anions), AGG-2 (three or four FSI⁻ anions), and AGG-3 (five or six FSI⁻ anions). In the mixed EC/DEC electrolyte, the solvation structures are classified as solvent-separated ion pairs (SSIPs), contact ion pairs (CIPs), and aggregates (AGGs). The solvation structures in DPE + 1.8 M LiFSI, FEME + 1.8 M LiFSI, and EC/DEC + 1 M LiPF₆ electrolytes are illustrated in Fig. S15–S19 (ESI†). The frequency of occurrence of all possible unique solvation structures in each electrolyte is analyzed in Fig. 9 and 10. Fig. 10b also shows the percentage of SSIP, CIP, and AGG in mixed EC/DEC electrolyte. The two most frequent solvation structures in each electrolyte and their respective percentages are shown in Table 4. In DPE + 1.8 M LiFSI electrolyte, the two most dominant solvation structures are Li⁺(FSI⁻)₃(DPE)₁ (48.15%) and Li⁺(FSI⁻)₄(DPE)₁ (17.78%). Similarly, in FEME + 1.8 M LiFSI electrolyte, the two most dominant solvation structures are Li⁺(FSI⁻)₃(FEME)₁ (33.33%) and Li⁺(FSI⁻)₄(FEME)₁ (17.78%). According to our study, the primary solvation structure remains nearly unchanged across varying salt concentrations in each DPE- and FEME-based electrolyte. Table 5 provides a thorough analysis of each solvation structure in the mixed EC/DEC electrolyte, specifying the individual counts of EC and DEC solvent molecules. Our findings also indicate that the number of unique solvation structures in DPE + LiFSI and FEME + LiFSI electrolytes increases as salt concentration increases.⁹⁵ In DPE + LiFSI, as the salt concentration rises from 1 M to 1.8 M to 4 M, the number of unique solvation structures increases from 7 to 8 to 10 (Fig. 9a–c). Similarly, in FEME + LiFSI, increasing the concentration from 1 M to 1.8 M to 4 M results in an increase in unique solvation structures from 9 to 10 to 12 (Fig. 9d–f). However, each FEME + LiFSI electrolyte exhibits a slightly higher number of possible solvation structures than DPE + LiFSI electrolytes. A detailed speciation of AGG in each DPE- and FEME-based electrolyte is shown in Fig. 11. At all salt concentrations, in DPE + LiFSI and FEME + LiFSI electrolytes, the primary aggregates belong to the AGG-2 category, containing three or four FSI⁻ anions and one solvent molecule (Fig. 11 and 12). Our results show that in each DPE- and FEME-based electrolyte, the solvation structures are primarily AGG (Fig. 9), whereas in the mixed EC/DEC electrolyte, 68% are SSIP and only 5.33% are AGG (Fig. 10b).

Fig. 9 Percentage of unique Li⁺ solvation structures in each fluorinated electrolyte: (a) DPE + 1 M LiFSI, (b) DPE + 1.8 M LiFSI, (c) DPE + 4 M LiFSI, (d) FEME + 1 M LiFSI, (e) FEME + 1.8 M LiFSI, and (f) FEME + 4 M LiFSI. The Li⁺ is surrounded by FSI⁻, DPE and FSI⁻, FEME.

Fig. 10 (a) Percentage of unique Li⁺ solvation structures in mixed carbonate 1 : 1 EC/DEC + 1 M LiPF₆ electrolyte. The Li⁺ is surrounded by EC, DEC, PF₆⁻. A value of PF₆ = 0 indicates that the PF₆⁻ anion is located outside the primary solvation shell (SSIPs). Inside the primary solvation shell, PF₆ = 1 refers to CIPs, and PF₆ > 1 refers to AGGs. (b) Percentage of SSIPs, CIPs, and AGGs.

Table 4 Top two most dominant solvation structures and their occurrence rates, calculated from MD simulations, in fluorinated and carbonate-based electrolytes with LiFSI and LiPF₆ salt at 25 °C

Electrolyte system	Solvation structure	Percentage (%)
DPE + 1 M LiFSI	Li^+(FSI^−)3(DPE)1	56.00
	Li^+(FSI^−)5(DPE)0	14.67
DPE + 1.8 M LiFSI	Li^+(FSI^−)3(DPE)1	48.15
	Li^+(FSI^−)4(DPE)1	17.78
DPE + 4 M LiFSI	Li^+(FSI^−)3(DPE)1	47.18
	Li^+(FSI^−)4(DPE)0	17.61
FEME + 1 M LiFSI	Li^+(FSI^−)5(FEME)0	24.00
	Li^+(FSI^−)3(FEME)1	22.67
FEME + 1.8 M LiFSI	Li^+(FSI^−)3(FEME)1	33.33
	Li^+(FSI^−)4(FEME)1	17.78
FEME + 4 M LiFSI	Li^+(FSI^−)3(FEME)1	30.23
	Li^+(FSI^−)4(FEME)1	18.60
1 : 1 EC/DEC + 1 M LiPF6	Li^+(PF6^−)0(EC/DEC)5	56.00
	Li^+(PF6^−)1(EC/DEC)4	16.00

---

Table 5 Solvation structures and their occurrence rates, calculated from MD simulations. A value of PF₆ = 0 indicates that the PF₆⁻ anion is located outside the primary solvation shell (SSIPs). Inside the primary solvation shell, PF₆ = 1 refers to CIPs, and PF₆ > 1 refers to AGGs

Solvation structures EC, DEC, PF6	Frequency of occurrence	Percentage (%)
(1, 4, 0)	25	33.33
(0, 5, 0)	10	13.33
(1, 3, 1)	7	9.33
(2, 3, 0)	7	9.33
(0, 3, 1)	5	6.67
(0, 4, 1)	4	5.33
(0, 4, 0)	3	4.00
(1, 3, 0)	3	4.00
(0, 2, 2)	3	4.00
(1, 2, 1)	2	2.67
(4, 2, 0)	2	2.67
(3, 3, 0)	1	1.33
(3, 2, 1)	1	1.33
(2, 2, 1)	1	1.33
(0, 3, 2)	1	1.33

---

Fig. 11 Detailed speciation of AGG in each fluorinated electrolyte with the frequency of occurrence. AGG-1 (one and two FSI⁻ anions), AGG-2 (three and four FSI⁻ anions), and AGG-3 (five and six FSI⁻ anions).

Fig. 12 Frequency of occurrence of the Li⁺ solvation structures in each fluorinated electrolyte based on the presence of zero, one, and two solvent molecules.

Li⁺ diffusion in electrolytes mainly happens through two types of mechanisms: the vehicle mechanism and the hopping mechanism.^121,122 In the vehicle mechanism, which is common in electrolytes with moderate salt concentration, Li⁺ ions move through the electrolyte where the Li⁺ ion and its coordinated solvent molecules migrate together as a complex. This can usually leads to a lower Li⁺ transference number because the ions movement is tied to the motion of the solvent and anions. In contrast, hopping mechanism is prevalent in concentrated electrolytes.^123–128 Here, Li⁺ ions move by jumping between different coordination sites in the electrolyte. These sites can be either solvent molecules or anions. This hopping mechanism often leads to higher Li⁺ mobility and transference numbers, as Li⁺ ions move more independently from the bulk solvent. Our MD and DFT studies show clear differences in Li⁺ diffusion mechanism depending on the solvation structure. The Li⁺ ions in the AGG-dominated DPE and FEME systems are likely to favor a hopping-type diffusion mechanism.^121,122 According to Saito et al., the increased microviscosity in WSEE, caused by strong cation–anion and cation–polymer interactions, slows down ion diffusion.¹²⁹ This, along with aggregate formation via bridging coordination, leads to suppressed ionic transport and reduced conductivity in the DPE electrolyte as confirmed by Li et al.⁴⁵ Conversely, in the SSIP-dominated mixed EC/DEC system, the Li⁺ ions are strongly solvated by solvent molecules and likely diffuse via a vehicle mechanism.^121,122 These differences in solvation structure and diffusion mechanism directly influence ionic conductivity, Li⁺ transference numbers, and ultimately battery performance. In future work, more electrolyte properties such as ionic conductivity, voltage window, and Li⁺ transference number can be explored.^45,130 These properties will further validate the transport behavior suggested by the observed solvation structures.

The charge density difference, Bader charge analysis, electrostatic potential maps, and binding energies of the solvation structures are discussed in the following sections.

3.3.1 Charge density difference and Bader charge analysis with ESP maps. In addition to Bader atomic charge calculations in VASP, the charge density difference (CDD) was also calculated for all the Li⁺ solvation structures using eqn (10).¹³¹ Here, ρ_totalsystem, ρ_subsystem1, and ρ_subsystem2 represent the charge densities of the Li⁺ solvation structure (Li + solvents + anions), Li, and the combined solvents and anions, respectively. These calculations, including CDD and Bader charge analysis, were performed to compute the amount of charge transfer from the Li⁺ ion to the surrounding solvents and anions and to analyze their electronic distribution.⁸⁸Fig. 13 shows the CDD plot of Li⁺ solvation structures, where the yellow regions indicate electron gain and the cyan regions represent electron loss. These CDD plots depict the regions of electron loss around the Li⁺ ion and electron accumulation between the ionized Li atom and the O atoms, indicating charge transfer from the ionized Li atom to the O atoms. This suggests a strong binding interaction between Li and the O atoms of the anions and solvents in the Li⁺ solvation structures. In the case of CIPs and AGGs in the EC/DEC electrolyte, strong binding interactions occur between Li and both the F atom of the PF₆⁻ anion and the O atom of the solvents. The amount of charge transfer from the Li⁺ ion to the surrounding solvents and anions was calculated using Bader charge analysis and is listed in Tables 6–10. Reduced charge transfer is observed in Li⁺–carbonate systems compared to Li⁺–ether systems (Table 6). Approximately 0.70–0.86e (Table 7) and 0.67–0.79e (Table 8) are transferred from Li⁺ to the surrounding solvents and anions in the AGGs of the 1.8 M DPE and 1.8 M FEME electrolytes. Similarly, approximately 0.79–0.85e, 0.71–0.84e, and 0.80e (Table 9) are transferred from Li⁺ to the surrounding solvents and anions in the SSIPs, CIPs, and AGGs of the 1 M EC/DEC electrolyte. Additionally, the ESP maps of these Li⁺ solvation structures in Fig. 13 show that the negative charge is primarily localized on the O and F atoms, while the positive charge is mainly localized on the Li atom.

Δρ = ρ_totalsystem − ρ_subsystem1 − ρ_subsystem2 (10)

Fig. 13 Top two most solvated structures (Li⁺(FSI⁻)₃(DPE)₁, Li⁺(FSI⁻)₄(DPE)₁, Li⁺(FSI⁻)₃(FEME)₁, Li⁺(FSI⁻)₄(FEME)₁, Li⁺(PF₆⁻)₀(EC)₁(DEC)₄, and Li⁺(PF₆⁻)₀(EC)₀(DEC)₅ including Li⁺(FSI⁻)₄(DPE)₀, Li⁺(FSI⁻)₄(FEME)₀), CIP, and AGG with their HOMO/LUMO and corresponding ESP and CDD in (a), (d) DPE + 1.8 M LiFSI, (b) and (e) FEME + 1.8 M LiFSI, and (c) and (f) EC/DEC + 1 M LiPF₆ electrolytes. The positive and negative phase of HOMO and LUMO are depicted in yellow and cyan colors, respectively. Yellow and cyan indicate the different signs of the isosurface of the wave function, and their sizes indicate its amplitude. In the ESP maps, the red and blue regions represent areas of high electron density (negative charge) and low electron density (positive charge), respectively. In the CDD plot, the cyan region represents electron depletion and the yellow region represents electron accumulation. The isosurfaces of HOMO/LUMO, ESP, and CDD were visualized using VESTA. The isosurface levels were set between 1 × 10⁻¹⁰ and 1 × 10⁻⁸ for HOMO/LUMO, 10 for ESP, and 0.0009 for CDD. HOMO/LUMO diagrams are shown in (a)–(c), while ESP maps and CDD plots are presented in (d)–(f).

Table 6 Calculated quantities of Li⁺–solvent and Li⁺–anion systems

Systems	HOMO (eV)	LUMO (eV)	LUMO–HOMO (eV)	Binding energy (eV)	Bader charge of Li, q(e)	Δq(e)
Li^+(DPE)1	−1.687	−0.638	1.05	−0.561	+0.726	+0.274
Li^+(FEME)1	−2.002	−0.806	1.20	−0.531	+0.730	+0.270
Li^+(FSI^−)1	−6.270	−0.993	5.28	−5.431	+0.160	+0.840
Li^+(EC)1	−1.742	−1.394	0.35	−0.651	+0.788	+0.212
Li^+(DEC)1	−1.698	−1.156	0.54	−0.671	+0.801	+0.199
Li^+(PF6^−)1	−8.077	−1.158	6.92	−7.491	+0.053	+0.947

---

Table 7 Calculated quantities of Li⁺ solvation structures in DPE + 1.8 M LiFSI electrolyte

Solvation structure	Species	HOMO (eV)	LUMO (eV)	LUMO–HOMO (eV)	Binding energy (eV)	Bader charge of Li, q(e)	Δq(e)
Li^+(FSI^−)2(DPE)1	AGG1	−6.949	−0.940	6.01	−7.901	+0.146	+0.854
Li^+(FSI^−)2(DPE)2	AGG1	−6.076	−0.482	5.59	−8.331	+0.141	+0.859
Li^+(FSI^−)3(DPE)0	AGG2	−7.279	−1.554	5.73	−7.291	+0.202	+0.798
Li^+(FSI^−)4(DPE)0	AGG2	−6.981	−1.350	5.63	−7.501	+0.195	+0.805
Li^+(FSI^−)3(DPE)1	AGG2	−6.590	−0.966	5.62	−8.011	+0.299	+0.701
Li^+(FSI^−)4(DPE)1	AGG2	−6.328	−0.838	5.49	−8.601	+0.235	+0.765
Li^+(FSI^−)3(DPE)2	AGG2	−6.118	−0.406	5.71	−8.871	+0.286	+0.714
Li^+(FSI^−)5(DPE)0	AGG3	−6.364	−2.702	3.66	−8.911	+0.251	+0.749

---

Table 8 Calculated quantities of Li⁺ solvation structures in FEME + 1.8 M LiFSI electrolyte

Solvation structure	Species	HOMO (eV)	LUMO (eV)	LUMO–HOMO (eV)	Binding energy (eV)	Bader charge of Li, q (e)	Δq(e)
Li^+(FSI^−)2(FEME)1	AGG1	−6.927	−1.101	5.83	−7.621	+0.239	+0.761
Li^+(FSI^−)2(FEME)2	AGG1	−6.405	−0.875	5.53	−8.001	+0.232	+0.768
Li^+(FSI^−)3(FEME)0	AGG2	−7.200	−1.543	5.66	−7.361	+0.226	+0.774
Li^+(FSI^−)4(FEME)0	AGG2	−6.841	−1.365	5.48	−7.771	+0.309	+0.691
Li^+(FSI^−)3(FEME)1	AGG2	−6.787	−0.972	5.81	−7.761	+0.314	+0.686
Li^+(FSI^−)4(FEME)1	AGG2	−6.570	−0.915	5.66	−8.171	+0.252	+0.748
Li^+(FSI^−)3(FEME)2	AGG2	−6.448	−0.628	5.82	−8.381	+0.227	+0.773
Li^+(FSI^−)5(FEME)0	AGG3	−6.664	−1.037	5.63	−7.951	+0.266	+0.734
Li^+(FSI^−)6(FEME)0	AGG3	−6.406	−0.829	5.58	−8.261	+0.215	+0.785
Li^+(FSI^−)5(FEME)1	AGG3	−6.347	−0.719	5.63	−8.641	+0.330	+0.670

---

Table 9 Calculated quantities of Li⁺ solvation structures in EC/DEC + 1 M LiPF₆ electrolyte

Solvation structure	Species	HOMO (eV)	LUMO (eV)	LUMO–HOMO (eV)	Binding energy (eV)	Bader charge of Li, q(e)	Δq(e)
Li^+(EC)0(DEC)4	SSIP	−6.433	−0.860	5.57	−10.151	+0.214	+0.786
Li^+(EC)0(DEC)5	SSIP	−6.241	−0.808	5.43	−10.591	+0.158	+0.842
Li^+(EC)1(DEC)3	SSIP	−6.021	−1.216	4.81	−9.481	+0.192	+0.808
Li^+(EC)1(DEC)4	SSIP	−6.484	−0.667	5.82	−10.761	+0.192	+0.808
Li^+(EC)2(DEC)3	SSIP	−6.559	−0.786	5.77	−10.621	+0.203	+0.797
Li^+(EC)3(DEC)3	SSIP	−6.118	−0.856	5.26	−10.741	+0.149	+0.851
Li^+(EC)4(DEC)2	SSIP	−5.704	−1.076	4.63	−10.121	+0.164	+0.836
Li^+(PF6^−)1(EC)0(DEC)3	CIP	−6.988	−0.877	6.11	−9.871	+0.165	+0.835
Li^+(PF6^−)1(EC)0(DEC)4	CIP	−6.061	−0.736	5.33	−10.121	+0.289	+0.711
Li^+(PF6^−)1(EC)1(DEC)2	CIP	−6.765	−1.003	5.76	−9.751	+0.235	+0.765
Li^+(PF6^−)1(EC)1(DEC)3	CIP	−6.648	−0.759	5.89	−10.201	+0.160	+0.840
Li^+(PF6^−)1(EC)2(DEC)2	CIP	−6.459	−0.743	5.72	−10.201	+0.231	+0.769
Li^+(PF6^−)1(EC)3(DEC)2	CIP	−6.517	−0.719	5.80	−11.051	+0.177	+0.823
Li^+(PF6^−)2(EC)0(DEC)2	AGG	−7.808	−1.465	6.34	−10.971	+0.198	+0.802
Li^+(PF6^−)2(EC)0(DEC)3	AGG	−7.390	−1.104	6.29	−11.551	+0.196	+0.804

---

Table 10 Calculated quantities of SSIPs in EC/DEC + 1 M LiPF₆ electrolyte. PF₆⁻ anion is removed from SSIPs

Solvation structure	Species	HOMO (eV)	LUMO (eV)	LUMO–HOMO (eV)	Binding energy (eV)	Bader charge of Li, q(e)	Δq(e)
Li^+(EC)0(DEC)4	SSIP	−0.912	−0.649	0.26	−2.611	+0.179	+0.821
Li^+(EC)0(DEC)5	SSIP	−0.756	−0.565	0.19	−3.061	+0.151	+0.849
Li^+(EC)1(DEC)3	SSIP	−0.947	−0.767	0.18	−2.391	+0.207	+0.793
Li^+(EC)1(DEC)4	SSIP	−0.867	−0.599	0.27	−3.001	+0.187	+0.813
Li^+(EC)2(DEC)3	SSIP	−0.891	−0.643	0.25	−2.891	+0.189	+0.811
Li^+(EC)3(DEC)3	SSIP	−0.838	−0.536	0.30	−3.461	+0.176	+0.824
Li^+(EC)4(DEC)2	SSIP	−0.793	−0.526	0.27	−3.081	+0.204	+0.796

---

3.3.2 Binding energy. The binding energies of different lithium-ion solvated systems were investigated using DFT calculations. The binding energy of the Li⁺ solvation structure was calculated by subtracting the energies of individual components from the total energy of the Li⁺ solvation structure, using eqn (11) and (12), where E represents the energy, and n and m are the number of solvent and anion species in the lithium-ion solvation structure, respectively.^45,95 Our findings show that DPE and FEME exhibit significantly reduced binding energy to Li⁺ (Fig. 14 and Table 6), which effectively suppresses salt dissociation and promotes ion aggregate formation starting at a concentration of 1 M.⁴⁵ The binding energies of these AGGs in the 1.8 M DPE and 1.8 M FEME electrolytes are quite similar, ranging from −8.91 to −7.29 eV and −8.64 to −7.36 eV, respectively (Tables 7 and 8). For the SSIPs (including PF₆⁻), CIPs, and AGGs, the binding energy ranges from −11.55 to −9.48 eV, while for SSIPs (excluding PF₆⁻), the binding energy is lower, ranging from −3.46 to −2.39 eV (Tables 9 and 10). The DFT validation is shown in Fig. 14. The slight differences in binding energies for the same systems depend on the software used (VASP or Gaussian), the functional (PBE or B3LYP), the simulation box size, and the orientation of the structure.^{45,49,50,132–137}

Fig. 14 DFT validation of binding energies.^{45,49,50,132–137}

Binding energy is the energy required to break a system into its individual components and separate them infinitely. In the context of lithium-ion solvation structures, binding energy refers to the strength of interaction between the Li⁺ ion and the solvent or anions species in the cation solvation shell. It quantifies how tightly the Li⁺ ion is bound to its surrounding environment, including solvent molecules (DPE, FEME, EC, DEC) and the counterions (FSI⁻, PF₆⁻). If the binding energy is negative, bond formation is likely exergonic. Moreover, a higher absolute value of the binding energy reflects stronger interactions among the species.⁹⁴ According to the literature, during the desolvation process, Li⁺ ions separate from the solvated molecules, move through the solid electrolyte interphase (SEI) layer, and diffuse into the electrode.¹³⁸ The binding energy in the Li⁺ solvation structure is important because it affects how easily the Li⁺ ion moves through the electrolyte and intercalates into the electrode, which influences battery performance. Recent studies have shown that using isobutyronitrile (iBN) as a cosolvent weakens the Li⁺–solvent interaction, making Li⁺ desolvation easier and thereby improving low-temperature ionic mobility.¹³⁹ A higher binding energy (greater stability of solvation) means stronger interactions between the Li⁺ ion and surrounding molecules, making it harder to desolvate the Li⁺ ion and slowing battery performance, especially at low temperatures.¹⁴⁰ Conversely, lower binding energy allows the Li⁺ ion to move more easily, improving battery efficiency. From Fig. 14 and Table 6, the binding energy values indicate that DEC (−0.671 eV) and EC (−0.651 eV) exhibit stronger solvating power than DPE (−0.561 eV) and FEME (−0.531 eV), as more negative binding energies correspond to stronger interactions with Li⁺. This trend is also consistent with the higher coordination with Li⁺ observed in carbonate solvents compared to ether solvents (Table 3). Furthermore, among the ethers, FEME exhibits even weaker solvating power than DPE, as reflected by its less negative binding energy.

E_b = E_complex − (E_Li⁺ + nE_solvent + mE_anion) (11)

E_b = E_complex − (E_Li⁺ + n₁E_solvent1 + n₂E_solvent2 + mE_anion) (12)

3.3.3 HOMO and LUMO distributions. The HOMO and LUMO distributions of all the lithium-ion solvation structures were simulated to investigate their reductive stability and to understand their decomposition mechanisms (Fig. 13). The LUMO energy level often serves as a key indicator of the reductive stability of electrolyte solutions. Additionally, the molecular orbital diagram of the LUMO can help determine potential decomposition pathways of these solvation structures, as the LUMO is the orbital where electron acceptance occurs during reduction.¹⁰⁰ Furthermore, the energy band gap (LUMO−HOMO) can also determine the chemical reactivity and stability of these solvation structures. A smaller band gap typically corresponds to higher chemical reactivity and lower stability, while a larger band gap suggests reduced reactivity and increased stability.⁵⁹ Among these solvation structures, the SSIPs without the PF₆⁻ anion have a lower band gap (0.18 to 0.30 eV) and binding energy (−3.50 to −2.40 eV), making them the most unstable structures (Table 10).

The LUMO is primarily distributed on the FSI⁻ anion in all LiFSI containing solvation structures in both 1.8 M DPE and 1.8 M FEME electrolytes (Fig. 13 and Fig. S20, S21, ESI†). Hence, the FSI⁻ anions will preferentially undergo reductive decomposition.¹⁰⁰ In the SSIP, CIP, and AGG structures in EC/DEC + 1 M LiPF₆ electrolyte, the LUMO is distributed across both EC and DEC molecules. These findings indicate that reductive decomposition reactions may occur through both EC and DEC decomposition.¹⁰⁰

4. Conclusions

In the present study, we systematically investigate the electrolyte structures in fluorinated ether (FEME + LiFSI), non-fluorinated ether (DPE + LiFSI), and organic carbonate-based (EC/DEC + LiPF₆) electrolytes over a wide range of salt concentrations (1 M, 1.8 M, and 4 M) using a combination of classical MD simulations with the OPLS-AA force field and DFT calculations. We observe that AGGs are the predominant species in the ether-based electrolytes, whereas SSIPs dominate in the mixed carbonate-based electrolyte. This aggregation effect is particularly strong in FEME-based electrolytes, supported by the high coordination number of Li⁺–Li⁺ pairs and the comparatively lower binding energy of FEME to Li⁺. The most dominant solvation structure in each ether-based electrolyte is the anion-rich solvation structure Li⁺(FSI⁻)₃(DPE)₁ and Li⁺(FSI⁻)₃(FEME)₁, respectively, and remain nearly unchanged across varying salt concentrations. Regarding the solvent composition in the solvation structures of the EC/DEC electrolyte, a higher fraction of DEC appears to be favorable. Our findings indicate that both DPE and FEME solvents exhibit weak solvating power at all salt concentrations, as indicated by the radial distribution functions, coordination numbers, and solvation structures, which show a strong preference for Li⁺ to interact with FSI⁻ anions in the primary solvation shell. In particular, FEME shows even weaker solvating power than DPE, as indicated by the higher coordination numbers of FSI⁻ in the primary solvation shell of FEME electrolytes. We also observe an increase in unique solvation structures in ether-based electrolytes with higher salt concentrations, with FEME + LiFSI displaying a slightly larger variety of structures than DPE + LiFSI. Furthermore, the electronic information of the lithium-ion solvation structures obtained from the DFT calculations are quite similar for both DPE- and FEME-based electrolytes. The charge density difference and Bader charge analysis show that the charge transfer from Li⁺ to the surrounding solvents and anions in the AGGs of the DPE + 1.8 M LiFSI electrolyte (0.70–0.86e) is comparatively higher than in the FEME + 1.8 M LiFSI electrolyte (0.67–0.79e). The binding energies of these AGGs in the 1.8 M DPE and 1.8 M FEME electrolytes are quite similar, ranging from −8.91 to −7.29 eV and −8.64 to −7.36 eV, respectively. The chemical stability of the solvation structures has also been predicted using their HOMO/LUMO distributions. Fluorinated electrolytes present safety concerns, including volatility and flammability.¹⁴¹ This study computationally investigates the solvation structure of FEME electrolyte and does not propose for immediate commercial use due to the risk of flammability. Further experimental investigation is needed to gain a deeper understanding of anion-rich solvation structures and flammability in FEME-based electrolytes. Combined with experimental and computational studies, our findings could provide valuable insights for advancing AGG-dominated FEME-based electrolyte design to meet the demands of LiF-rich SEI layers in next-generation lithium-ion batteries.

Author contributions

Conceptualization: R. H.; simulation (MD and DFT): R. H.; data management: R. H.; compilation: R. H.; proofreading: R. H. and D. D.; supervision: D. D.; funding acquisition: D. D.

Conflicts of interest

There are no conflicts to declare.

Data availability

The input scripts and data files for running MD simulations in LAMMPS, along with the optimized structure files (DPE, FEME, EC, DEC, LiFSI, and LiPF₆) from VASP, are available on the GitHub repository at https://github.com/mana121/SolvationStructure.git. The ESI† includes GIF generated using the OVITO software, illustrating the full 5 ns trajectory of the production run for the DPE + 1.8 M LiFSI, FEME + 1.8 M LiFSI, and 1 : 1 EC/DEC + 1 M LiPF₆ electrolyte systems.

Acknowledgements

This work used resources of the HPC cluster Wulver at the New Jersey Institute of Technology (NJIT) and Expanse, a dedicated Advanced Cyberinfrastructure Coordination Ecosystem: Services and Support (ACCESS) cluster, which are supported by National Science Foundation (NSF) award numbers CMMI-2237990, CBET-2126180, and ACCESS awarded resources or project DMR180013. The authors extend their appreciation to the ViPER Lab of Vilas Pol at Purdue University for their valuable observations on the initial simulation results of electrolytes conducted by the author.

References

1. K. Xu, Chem. Rev., 2004, 104, 4303–4418.
2. M. Armand and J.-M. Tarascon, Nature, 2008, 451, 652–657.
3. J. B. Goodenough and Y. Kim, Chem. Mater., 2010, 22, 587–603.
4. J.-M. Tarascon and M. Armand, Nature, 2001, 414, 359–367.
5. B. Scrosati, J. Hassoun and Y.-K. Sun, Energy Environ. Sci., 2011, 4, 3287–3295.
6. M. Li, J. Lu, Z. Chen and K. Amine, Adv. Mater., 2018, 30, 1800561.
7. Y. Wu, L. Xie, H. Ming, Y. Guo, J.-Y. Hwang, W. Wang, X. He, L. Wang, H. N. Alshareef and Y.-K. Sun, et al. , ACS Energy Lett., 2020, 5, 807–816.
8. J. Xu, J. Zhang, T. P. Pollard, Q. Li, S. Tan, S. Hou, H. Wan, F. Chen, H. He and E. Hu, et al. , Nature, 2023, 614, 694–700.
9. L. Cheng, R. S. Assary, X. Qu, A. Jain, S. P. Ong, N. N. Rajput, K. Persson and L. A. Curtiss, J. Phys. Chem. Lett., 2015, 6, 283–291.
10. X. Qu, A. Jain, N. N. Rajput, L. Cheng, Y. Zhang, S. P. Ong, M. Brafman, E. Maginn, L. A. Curtiss and K. A. Persson, Comput. Mater. Sci., 2015, 103, 56–67.
11. O. Borodin, M. Olguin, C. E. Spear, K. W. Leiter and J. Knap, Nanotechnology, 2015, 26, 354003.
12. G. Bieker, M. Winter and P. Bieker, Phys. Chem. Chem. Phys., 2015, 17, 8670–8679.
13. Y.-C. Tsai, C.-T. Kuo, S.-F. Liu, Y.-T. Lee and T.-R. Yew, J. Phys. Chem. C, 2021, 125, 1221–1233.
14. M. Baek, J. Kim, J. Jin and J. W. Choi, Nat. Commun., 2021, 12, 6807.
15. Y. Chae, C. Lim, J. Jeon, M. Kim, K.-K. Lee, K. Kwak and M. Cho, J. Phys. Chem. Lett., 2022, 13, 7881–7888.
16. J. Lim, K.-K. Lee, C. Liang, K.-H. Park, M. Kim, K. Kwak and M. Cho, J. Phys. Chem. B, 2019, 123, 6651–6663.
17. K. Fulfer and D. Kuroda, Phys. Chem. Chem. Phys., 2017, 19, 25140–25150.
18. K.-K. Lee, K. Park, H. Lee, Y. Noh, D. Kossowska, K. Kwak and M. Cho, Nat. Commun., 2017, 8, 14658.
19. K. D. Fulfer and D. G. Kuroda, J. Phys. Chem. C, 2016, 120, 24011–24022.
20. D. M. Seo, S. Reininger, M. Kutcher, K. Redmond, W. B. Euler and B. L. Lucht, J. Phys. Chem. C, 2015, 119, 14038–14046.
21. N. Chapman, O. Borodin, T. Yoon, C. C. Nguyen and B. L. Lucht, J. Phys. Chem. C, 2017, 121, 2135–2148.
22. A. V. Cresce, S. M. Russell, O. Borodin, J. A. Allen, M. A. Schroeder, M. Dai, J. Peng, M. P. Gobet, S. G. Greenbaum and R. E. Rogers, et al. , Phys. Chem. Chem. Phys., 2017, 19, 574–586.
23. S. Hwang, D.-H. Kim, J. H. Shin, J. E. Jang, K. H. Ahn, C. Lee and H. Lee, J. Phys. Chem. C, 2018, 122, 19438–19446.
24. H. Lee, S. Hwang, M. Kim, K. Kwak, J. Lee, Y.-K. Han and H. Lee, J. Phys. Chem. Lett., 2020, 11, 10382–10387.
25. X. Bogle, R. Vazquez, S. Greenbaum, A. V. W. Cresce and K. Xu, J. Phys. Chem. Lett., 2013,(4), 1664–1668.
26. I. Skarmoutsos, V. Ponnuchamy, V. Vetere and S. Mossa, J. Phys. Chem. C, 2015, 119, 4502–4515.
27. O. Borodin, M. Olguin, P. Ganesh, P. R. Kent, J. L. Allen and W. A. Henderson, Phys. Chem. Chem. Phys., 2016, 18, 164–175.
28. Y.-Q. Wang, H. Xu, B. Cao, J. Ma and Z.-W. Yu, J. Phys. Chem. Lett., 2024, 15, 5047–5055.
29. P. Wróbel, P. Kubisiak and A. Eilmes, J. Phys. Chem. B, 2021, 125, 1248–1258.
30. O. Borodin and G. D. Smith, J. Phys. Chem. B, 2009, 113, 1763–1776.
31. C. Liang, K. Kwak and M. Cho, J. Phys. Chem. Lett., 2017, 8, 5779–5784.
32. S. Han, Sci. Rep., 2019, 9, 5555.
33. M. Mynam, B. Ravikumar and B. Rai, J. Mol. Liq., 2019, 278, 97–104.
34. X. Zhang and D. G. Kuroda, J. Chem. Phys., 2019, 150, 184501.
35. M. J. Boyer, L. Vilciauskas and G. S. Hwang, Phys. Chem. Chem. Phys., 2016, 18, 27868–27876.
36. T. Hou, K. D. Fong, J. Wang and K. A. Persson, Chem. Sci., 2021, 12, 14740–14751.
37. H. Cheng, Q. Sun, L. Li, Y. Zou, Y. Wang, T. Cai, F. Zhao, G. Liu, Z. Ma and W. Wahyudi, et al. , ACS Energy Lett., 2022, 7, 490–513.
38. Z. Piao, R. Gao, Y. Liu, G. Zhou and H.-M. Cheng, Adv. Mater., 2023, 35, 2206009.
39. Y.-X. Yao, X. Chen, C. Yan, X.-Q. Zhang, W.-L. Cai, J.-Q. Huang and Q. Zhang, Angew. Chem., Int. Ed., 2021, 60, 4090–4097.
40. T. D. Pham and K.-K. Lee, Small, 2021, 17, 2100133.
41. J. Chen, X. Fan, Q. Li, H. Yang, M. R. Khoshi, Y. Xu, S. Hwang, L. Chen, X. Ji and C. Yang, et al. , Nat. Energy, 2020, 5, 386–397.
42. K. Chen, X. Shen, L. Luo, H. Chen, R. Cao, X. Feng, W. Chen, Y. Fang and Y. Cao, Angew. Chem., 2023, 135, e202312373.
43. C.-C. Su, M. He, R. Amine, T. Rojas, L. Cheng, A. T. Ngo and K. Amine, Energy Environ. Sci., 2019, 12, 1249–1254.
44. C. Reichardt, J. Org. Chem., 2021, 87, 1616–1629.
45. Z. Li, H. Rao, R. Atwi, B. M. Sivakumar, B. Gwalani, S. Gray, K. S. Han, T. A. Everett, T. A. Ajantiwalay and V. Murugesan, et al. , Nat. Commun., 2023, 14, 868.
46. K. Xu, Chem. Rev., 2014, 114, 11503–11618.
47. K. Xu, A. von Cresce and U. Lee, Langmuir, 2010, 26, 11538–11543.
48. Y. Liu, Y. Zhu and Y. Cui, Nat. Energy, 2019, 4, 540–550.
49. Z. Wang, Z. Sun, Y. Shi, F. Qi, X. Gao, H. Yang, H.-M. Cheng and F. Li, Adv. Energy Mater., 2021, 11, 2100935.
50. W. Zhang, T. Yang, X. Liao, Y. Song and Y. Zhao, Energy Storage Mater., 2023, 57, 249–259.
51. Y. Zou, F. Cheng, Y. Lu, Y. Xu, C. Fang and J. Han, Small, 2023, 19, 2203394.
52. Y. Zhao, T. Zhou, M. Mensi, J. W. Choi and A. Coskun, Nat. Commun., 2023, 14, 299.
53. Z. Piao, H.-R. Ren, G. Lu, K. Jia, J. Tan, X. Wu, Z. Zhuang, Z. Han, C. Li and R. Gao, et al. , Angew. Chem., 2023, 135, e202300966.
54. T. Zhou, Y. Zhao, M. El Kazzi, J. W. Choi and A. Coskun, Angew. Chem., Int. Ed., 2022, 61, e202115884.
55. H. Chen, K. Chen, L. Luo, X. Liu, Z. Wang, A. Zhao, H. Li, X. Ai, Y. Fang and Y. Cao, Angew. Chem., Int. Ed., 2024, 63, e202316966.
56. Ø. Gullbrekken, A. F. Gunnarshaug, A. Lervik, S. Kjelstrup and S. K. Schnell, J. Am. Chem. Soc., 2024, 146, 4592–4604.
57. J. Wu, X. Wang, Q. Liu, S. Wang, D. Zhou, F. Kang, D. Shanmukaraj, M. Armand, T. Rojo and B. Li, et al. , Nat. Commun., 2021, 12, 5746.
58. Z. Jiang, T. Yang, C. Li, J. Zou, H. Yang, Q. Zhang and Y. Li, Adv. Funct. Mater., 2023, 33, 2306868.
59. X. He, Y. Zhang, K. Wang, B. Shan, M. Zhou, W. Wang and K. Jiang, J. Mater. Chem. A, 2022, 10, 23068–23078.
60. N. Vangapally, D. Lusztig, S. Rathod, A. Bano, H. Scalar, S. H. Akella, M. Noked, D. T. Major, I. C. Halalay and S. Sriramulu, et al. , Small, 2024, 20, 2403694.
61. W.-h Hou, Q. Feng, C. Liu, X. Zhang, J. Yue, Q. Tian, S. Wu, Y. Ou, P. Zhou and Y. Xia, et al. , Adv. Mater., 2025, 2503893.
62. Z. Jiang, J. Mo, C. Li, H. Li, Q. Zhang, Z. Zeng, J. Xie and Y. Li, Energy Environ. Mater., 2023, 6, e12440.
63. Q. Ma, S. Fu, A.-J. Wu, Q. Deng, W.-D. Li, D. Yue, B. Zhang, X.-W. Wu, Z.-L. Wang and Y.-G. Guo, Adv. Energy Mater., 2023, 13, 2203892.
64. Q. Wang, Z. Yao, C. Zhao, T. Verhallen, D. P. Tabor, M. Liu, F. Ooms, F. Kang, A. Aspuru-Guzik and Y.-S. Hu, et al. , Nat. Commun., 2020, 11, 4188.
65. P.-P. Chen, B.-H. Zhang, Z.-A. Li, J.-T. Lei, J.-Z. Chen, Y.-L. Hou and D.-L. Zhao, Small, 2024, 20, 2403079.
66. Z. Li, J. Liu, X. Bi, Y. Qin and T. Gao, J. Mater. Chem. A, 2023, 11, 19996–20010.
67. J. Moon, D. O. Kim, L. Bekaert, M. Song, J. Chung, D. Lee, A. Hubin and J. Lim, Nat. Commun., 2022, 13, 4538.
68. L. Martnez, R. Andrade, E. G. Birgin and J. M. Martnez, J. Comput. Chem., 2009, 30, 2157–2164.
69. W. L. Jorgensen, J. D. Madura and C. J. Swenson, J. Am. Chem. Soc., 1984, 106, 6638–6646.
70. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
71. W. L. Jorgensen and J. Tirado-Rives, Proc. Natl. Acad. Sci. U. S. A., 2005, 102, 6665–6670.
72. L. S. Dodda, J. Z. Vilseck, J. Tirado-Rives and W. L. Jorgensen, J. Phys. Chem. B, 2017, 121, 3864–3870.
73. L. S. Dodda, I. Cabeza de Vaca, J. Tirado-Rives and W. L. Jorgensen, Nucleic Acids Res., 2017, 45, W331–W336.
74. J. N. Canongia Lopes, J. Deschamps and A. A. Pádua, J. Phys. Chem. B, 2004, 108, 2038–2047.
75. K. Shimizu, D. Almantariotis, M. F. C. Gomes, A. A. Pádua and J. N. Canongia Lopes, J. Phys. Chem. B, 2010, 114, 3592–3600.
76. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In't Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al. , Comput. Phys. Commun., 2022, 271, 108171.
77. R. W. Hockney and J. W. Eastwood, Computer simulation using particles, CRC Press, 2021.
78. J.-P. Ryckaert, G. Ciccotti and H. J. Berendsen, J. Comput. Phys., 1977, 23, 327–341.
79. W. Shinoda, M. Shiga and M. Mikami, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 134103.
80. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695.
81. S. Nosé, Mol. Phys., 1984, 52, 255–268.
82. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
83. J. Harl, L. Schimka and G. Kresse, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 115126.
84. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
85. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
86. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
87. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
88. J. Zhang, Y.-F. Zhang, Y. Li, Y.-R. Ren, S. Huang, W. Lin and W.-K. Chen, Phys. Chem. Chem. Phys., 2021, 23, 5143–5151.
89. X. Kou, Y. Ma, C. Pan, Y. Huang, Y. Duan and Y. Yang, Langmuir, 2022, 38, 6116–6127.
90. S. Li, J. Zhang, S. Zhang, Q. Liu, H. Cheng, L. Fan, W. Zhang, X. Wang, Q. Wu and Y. Lu, Nat. Energy, 2024, 9, 285–297.
91. Y. Zhao, T. Zhou, T. Ashirov, M. E. Kazzi, C. Cancellieri, L. P. Jeurgens, J. W. Choi and A. Coskun, Nat. Commun., 2022, 13, 2575.
92. C. Wang, Z. Sun, Y. Liu, L. Liu, X. Yin, Q. Hou, J. Fan, J. Yan, R. Yuan and M. Zheng, et al. , Nat. Commun., 2024, 15, 6292.
93. M. Zhou, C. Ye and D. Xiang, Molecules, 2022, 27, 3266.
94. R. Atwi, M. Bliss, M. Makeev and N. N. Rajput, Sci. Rep., 2022, 12, 15760.
95. H. Yang, J. Ji, H. Li and B. Liang, Ind. Eng. Chem. Res., 2023, 62, 6516–6524.
96. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
97. V. Wang, N. Xu, J.-C. Liu, G. Tang and W.-T. Geng, Comput. Phys. Commun., 2021, 267, 108033.
98. Z. Yu, H. Wang, X. Kong, W. Huang, Y. Tsao, D. G. Mackanic, K. Wang, X. Wang, W. Huang and S. Choudhury, et al. , Nat. Energy, 2020, 5, 526–533.
99. Y. Wang, C. Zheng, W. Xie, X. Liu, Y. Lu, Y. Hou, T. Ma, Z. Yan and J. Chen, Adv. Mater., 2024, 36, 2312302.
100. Y. Aoki, M. Oda, S. Kojima, Y. Yamaga, T. Ishihama, T. Nagashima, T. Doi and M. Inaba, J. Phys. Chem. C, 2022, 127, 69–77.
101. Y. Wu, D. Ruan, Z. Cui, J. Fan, D. Wang and X. Ren, Chem. Sci., 2024, 15, 4238–4274.
102. L. Bai, B. Chen, D. Zhang, Y. Xu, S. Zhang, Y. He, F. Yu, K. Liu, D. Xu and Z. Chang, Energy Fuels, 2024, 38, 10634–10652.
103. H. Guo, Y. Tian, Y. Liu, Y. Bai, J. Wu, F. Kang and B. Li, ACS Appl. Mater. Interfaces, 2022, 15, 1201–1209.
104. G. Song, Z. Yi, F. Su, L. Xie, Z. Wang, X.-X. Wei, G. Xu and C.-M. Chen, ACS Energy Lett., 2023, 8, 1336–1343.
105. J. Guo, X. Liu, Z. Shen, Y. Lv, X. Zhang, C. Zhang and X. Zhang, Adv. Funct. Mater., 2024, 2405951.
106. D. Jin, C. Park, J. Han, C. B. Dzakpasu, E. Kim, J. Oh, K. M. Kim, S. Lee and Y. M. Lee, Battery Energy, 2023, 2, 20220034.
107. X. Yang, X. Liu, J. Han, Z. Liu and X. Zhang, J. Phys. Chem. B, 2023, 127, 3026–3040.
108. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2009, 18, 015012.
109. H. Lundgren, M. Behm and G. Lindbergh, J. Electrochem. Soc., 2014, 162, A413.
110. S.-I. Lee, U.-H. Jung, Y.-S. Kim, M.-H. Kim, D.-J. Ahn and H.-S. Chun, Korean J. Chem. Eng., 2002, 19, 638–644.
111. Y. R. Dougassa, C. Tessier, L. El Ouatani, M. Anouti and J. Jacquemin, J. Chem. Thermodyn., 2013, 61, 32–44.
112. C. A. Landry, J. G. Bergeron, O. Mathieu, K. L. Gering and J. C. Thomas, Batteries, 2025, 11, 44.
113. J. Holoubek, H. Liu, Z. Wu, Y. Yin, X. Xing, G. Cai, S. Yu, H. Zhou, T. A. Pascal and Z. Chen, et al. , Nat. Energy, 2021, 6, 303–313.
114. H. Li, Y. Kang, W. Wei, C. Yan, X. Ma, H. Chen, Y. Sang, H. Liu and S. Wang, Nano-Micro Lett., 2024, 16, 197.
115. X. Cao, P. Gao, X. Ren, L. Zou, M. H. Engelhard, B. E. Matthews, J. Hu, C. Niu, D. Liu and B. W. Arey, et al. , Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2020357118.
116. B. Koo, H. Lee, S. Hwang, J. Lee, Y.-K. Han, K. H. Ahn, C. Lee and H. Lee, J. Phys. Chem. C, 2023, 127, 18271–18278.
117. Y. Aoki, M. Oda, S. Kojima, T. Ishihama, T. Nagashima, T. Doi and M. Inaba, ACS Appl. Energy Mater., 2022, 5, 1085–1094.
118. H. Lee, H. An, H. Chang, M. Lee, S. Park, S. Lee, J. Kang, S. Byon, B. Koo and H. Lee, et al. , Energy Storage Mater., 2023, 63, 102995.
119. H. Jia, P. Gao, L. Zou, K. S. Han, M. H. Engelhard, Y. He, X. Zhang, W. Zhao, R. Yi and H. Wang, et al. , Chem. Mater., 2020, 32, 8956–8964.
120. J. L. Banks, H. S. Beard, Y. Cao, A. E. Cho, W. Damm, R. Farid, A. K. Felts, T. A. Halgren, D. T. Mainz and J. R. Maple, et al. , J. Comput. Chem., 2005, 26, 1752–1780.
121. P. Zhou, X. Zhang, Y. Xiang and K. Liu, Nano Res., 2023, 16, 8055–8071.
122. S. Krämer, D. Weintz, M. Winter, I. Cekic-Laskovic and M. Grünebaum, Encyclopedia, 2025, 5, 20.
123. K. Dokko, D. Watanabe, Y. Ugata, M. L. Thomas, S. Tsuzuki, W. Shinoda, K. Hashimoto, K. Ueno, Y. Umebayashi and M. Watanabe, J. Phys. Chem. B, 2018, 122, 10736–10745.
124. D. Yu, D. Troya, A. G. Korovich, J. E. Bostwick, R. H. Colby and L. A. Madsen, ACS Energy Lett., 2023, 8, 1944–1951.
125. Y. Ugata, S. Sasagawa, R. Tatara, K. Ueno, M. Watanabe and K. Dokko, J. Phys. Chem. B, 2021, 125, 6600–6608.
126. S. Kondou, M. L. Thomas, T. Mandai, K. Ueno, K. Dokko and M. Watanabe, Phys. Chem. Chem. Phys., 2019, 21, 5097–5105.
127. Y. Ugata, M. L. Thomas, T. Mandai, K. Ueno, K. Dokko and M. Watanabe, Phys. Chem. Chem. Phys., 2019, 21, 9759–9768.
128. S. R. Galle Kankanamge and D. G. Kuroda, J. Phys. Chem. B, 2020, 124, 1965–1977.
129. Y. Saito, S. Takeda, W. Morimura, R. Kuratani and S. Nishikawa, J. Phys. Chem. C, 2017, 121, 23926–23930.
130. M. Lorenz and M. Schönhoff, J. Phys. Chem. B, 2024, 128, 2782–2791.
131. I. Choudhuri and D. G. Truhlar, J. Chem. Theory Comput., 2020, 16, 5884–5892.
132. Y. Chen, Q. He, Y. Zhao, W. Zhou, P. Xiao, P. Gao, N. Tavajohi, J. Tu, B. Li and X. He, et al. , Nat. Commun., 2023, 14, 8326.
133. X. Chen, N. Yao, B.-S. Zeng and Q. Zhang, Fundam. Res., 2021, 1, 393–398.
134. S.-J. Yang, N. Yao, X.-Q. Xu, F.-N. Jiang, X. Chen, H. Liu, H. Yuan, J.-Q. Huang and X.-B. Cheng, J. Mater. Chem. A, 2021, 9, 19664–19668.
135. X. Chen, X.-Q. Zhang, H.-R. Li and Q. Zhang, Batteries Supercaps, 2019, 2, 128–131.
136. P. Zhou, Y. Ou, Q. Feng, Y. Xia, H. Zhou, W.-H. Hou, X. Song, Y. Lu, S. Yan and W. Zhang, et al. , Adv. Funct. Mater., 2025, 35, 2416800.
137. T. Chen, Z. Jin, Y. Liu, X. Zhang, H. Wu, M. Li, W. Feng, Q. Zhang and C. Wang, Angew. Chem., Int. Ed., 2022, 61, e202207645.
138. Z. Wang, H. Wang, S. Qi, D. Wu, J. Huang, X. Li, C. Wang and J. Ma, EcoMat, 2022, 4, e12200.
139. L. Luo, K. Chen, H. Chen, H. Li, R. Cao, X. Feng, W. Chen, Y. Fang and Y. Cao, Adv. Mater., 2024, 36, 2308881.
140. C. Fu, L. Xu, F. W. Aquino, A. V. Cresce, M. Gobet, S. G. Greenbaum, K. Xu, B. M. Wong and J. Guo, J. Phys. Chem. Lett., 2018, 9, 1739–1745.
141. J. Hou, L. Wang, X. Feng, J. Terada, L. Lu, S. Yamazaki, A. Su, Y. Kuwajima, Y. Chen and T. Hidaka, et al. , Energy Environ. Mater., 2023, 6, e12297.

---

Footnote

† Electronic supplementary information (ESI) available: See DOI: https://doi.org/10.1039/d5ya00154d

---