Theoretical insights into the dissolution of LiFSI in weakly and strongly solvating solvents

Abstract

Cyclopentyl methyl ether, CPME, a weakly solvating solvent for lithium bis(fluorosulfonyl)imide (LiFSI) electrolyte used in Li ion batteries, has been found experimentally to lead to the formation of a more stable solid electrolyte interface (SEI) compared to 1,2-dimethoxyethane (DME). In this study, force field based molecular dynamics simulations and density functional theory based ab initio calculations were used to study the interaction and clustering behaviour in CPME and DME based electrolytes at different concentrations. Treatment of the solvent polarization in force fields via the electronic continuum model reproduces results that largely agree with experimental data. Investigations of the solvent clustering around Li⁺ and the strength of the interaction between the anion and Li⁺ in solution sheds light on the concentration dependence of the formation of FSI⁻/Li⁺ aggregates, which are known to influence the stability of the SEI. In essence, the strength of Li⁺-DME interaction demands that Li⁺ first bind to all available DME molecules before binding to FSI⁻. Cluster survival time analysis further points to a more stable and dominant cluster Li⁺-(DME)₂ compared to the counter-ion pair FSI⁻/Li⁺ at low concentrations of LiFSI.

---

Introduction

Understanding the chemistry of lithium-ion batteries is essential for adoption of electric vehicles worldwide. For a wider adoption, there is a need to design batteries that can operate at wider temperature ranges in consideration of the climates of different markets. Current lithium-ion batteries are known to fail to hold charge at very low temperatures. In this regard, battery electrolytes play a significant role. Lithium-ion battery electrolytes are usually made of lithium salts dissolved in carbonate or ether-based solvents. Lithium hexafluorophosphate dissolved in ethylene carbonate is one of the most widely used electrolyte systems.¹ Lithium bis(fluorosulfonyl)imide (LiFSI) and lithium bis(trifluoromethanesulfonyl)imide (LiTFSI) have also been widely used as salts.² Ethylene carbonate (EC) is highly effective in dissociating lithium phosphate due to its high dielectric constant. In addition, the reduction of ethylene carbonate plays a role in the formation of the solid electrolyte interface (SEI).³ The SEI acts as a barrier to prevent electrolyte reduction which can lead to a short battery lifetime.⁴ A desirable SEI should be composed mainly of inorganic salts such as LiF which have high ionic conductivity but low electronic conductivity. A Li solvation sheath which is rich in solvent molecules obviously promotes an SEI which is rich in organics. A common strategy has been to increase the concentration of electrolytes solutions in solvents (high concentration electrolytes-HCE). This can lead to high viscosity electrolytes which is undesirable. It has been proposed that tuning the strength of interaction between the Li/solvent and anion/solvent is key in designing safe, efficient and long-lasting batteries. Thus, the nature of the solvent influences the stability of the solid electrolyte interface.⁵ Other solvents that have been used include dimethoxyethane (DME), 1,3-dioxolane (DOL), fluoroethylene carbonate (FEC), dimethyl carbonate (DMC), etc.² The interaction between DME and Li is primarily via the two oxygen atoms. For the monosolvated Li ion, ab initio calculations have shown that the binding energy between lithium and the solvent increases in the order FEC < EC < DOL < DMC < DEC < DME.⁶ Weaker Li–solvent interactions are generally preferred as desolvation step is usually the rate determining step in lithium-ion dynamics during charging.⁷

The solvating strength of ethers can be tuned in various ways including reducing ethereal groups, introduction of alkyl groups and steric hindrance.⁷ In a quest to find a solvent that can extend the liquid temperature range for lithium-ion batteries and in turn improve the coulombic efficiency, cyclopentyl methyl ether (CPME) was proposed.^8,9 CPME, an ether with a terminal methyl group, one ethereal oxygen and a 5-membered ring, is a weakly solvating solvent. Ramasamy et al. showed that Li/FSI aggregates are the dominant species in LiFSI/CPME electrolyte solutions even at 1 M LiFSI concentration unlike in LiFSI/DME electrolyte that is dominated by free FSI anions.⁸ The electrolyte had an initial efficiency of ∼87% and excellent low temperature performance. Similar work was reported by Zhang et al.⁹ showing that close to saturation (LiFSI to CPME ratio of 1 : 2) the electrolyte results in a high coulombic efficiency of >99% after 350 cycles on Li–Cu cell. In both papers, based on XPS data, it was shown that the SEI is dominated by inorganic salts in CPME based electrolytes. Recently, Wang and co-workers¹⁰ investigated how the performance of lithium–sulphur batteries was affected by the anodic stability of lithium metal and the cathodic conversion kinetics of sulphur–carbon composites. Their findings show that weakly solvating electrolytes like CPME, which moderately dissolves lithium polysulfides (LiPS), can suppress LiPS shuttling by forming an anion-dominated Li⁺ solvation structure. This effectively induces a stable solid electrolyte interphase (SEI), ensuring the anodic stability of the Li-metal anode. The behaviour is attributed to the low electronegativity of CPME and the steric hindrance introduced by its non-linear structure¹¹ unlike conventional DME which has high electron-donating ability and the ability to chelate.¹⁰ Replacing DME with the sterically hindered 1,2-dimethylpropane (DMP) as a solvent for 2 M LiFSI electrolytes also resulted in improved cycle stability of Li|LiFePO₄ and Li|LiNi_0.8Co_0.1Mn_0.1O₂ cells.¹² Force field based molecular dynamics simulations also traced the better performance of DMP based cells to the formation of Li/FSI aggregates. In this paper we apply systematic molecular dynamics simulation to study the LiFSI/CPME and LiFSI/DME electrolyte systems at different concentrations.

Computational methods

Molecular dynamics studies were conducted to investigate the solvation of lithium bis(fluorosulfonyl)imide (LiFSI) using solvents 1,2-dimethoxy ethane (DME) and cyclopentyl methyl ether (CPME) (see SI Fig. S1 for the structures). The construction of the configurations, the molecular dynamics simulations along with some of the analysis were performed using the GROMACS software.^13,14 The all-atom optimized potentials for liquids simulation (OPLS-AA) force field^15,16 was used for solvent molecules. The force field parameters for these solvents, including the charges, were obtained from the Ligpargen engine.¹⁷ The Canongia Lopes and Pádua (CL&P) force field for ionic liquids was used for the FSI anion.^18–20 The CL&P force field has the same functional form as the OPLS field. The Lennard-Jones parameters for Li ion are those used by Dang.²¹ For Li and FSI ions, the partial charges are scaled based on the Molecular Dynamics Electronic Continuum (MDEC) model.^22–25 The scaling accounts for electronic polarizability for solvents with low dielectric constants. For all partial charges in an ion, the effective partial charge, q^eff_i, is calculated as:²⁶where n is the refractive index of the solvent and q_i is the standard charge. Note that this scaling is applied to the charge of Li ion (+1) and all the partial charges in FSI⁻ only. The standard charges, q_i, in FSI⁻ are those from the CL&P force field. Initial configurations of the system were constructed by randomly inserting specified numbers of Li ions, bis(fluorosulfonyl)imide (FSI⁻) and solvents DME or CPME in a cubic box with dimensions of approximately 10 nm. The resulting mixtures were then energy minimized using the steepest descent algorithm. This was then followed by a 2 ns annealing under the following conditions. The energy minimized systems were gradually heated from room temperature to 348.15 K in 200 ps, where it was maintained for 200 ps. The temperature was then increased by 100 K to reach 448.15 K in 200 ps, and it was held there for 200 ps. Next, the temperature was raised by an additional 50 K to 498.15 K in 200 ps and maintained there for 200 ps. Finally, the system was cooled down to 293.15 K for the remaining isothermal isochoric dynamics (NVT). The systems were then subjected to isothermal isobaric dynamics (NPT) with T and P set to 298.15 K and 1 bar respectively, for 10–30 ns. The pressure was controlled through a Parrinello–Rahman barostat with compressibility set to 4.5 × 10⁻⁵ bar⁻¹ and a coupling constant 2 ps.²⁶ The temperature was controlled with a modified Berendsen thermostat²⁷ with a temperature coupling constant of 0.1 ps. The timestep for the simulations was set to 1 fs. The Verlet cut-off scheme^28,29 was used for the van der Waals and short-range electrostatic interactions and was kept at a distance of 1.2 nm. Long range electrostatic interactions were computed using the Particle Mesh Ewald algorithm.^30,31 Ab initio calculations were performed using Orca 6.1.0 quantum chemistry package.^32–35 Geometry optimizations and energy calculations were performed at the M06/6-31+G(d,p) level of theory. The MO6-2X version of the functional with 54% Hartree–Fock exchange was used.³⁶

Results and discussion

The interactions between solvents/anion/Li species in the electrolytes were studied via radial distribution functions (RDFs).

Shown in Fig. 1 are the RDFs between the FSI anion center of mass (COM) and the Li⁺ ion in different electrolytes at different concentrations. For CPME based electrolytes, the RDF peaks occur at approximately the same distance for different concentrations. However, the peak heights are different as expected for different concentrations. This indicates that the interaction mode between Li⁺ and FSI⁻ is the same at all three concentrations. Analysis of the three peaks below 500 pm shows that these are all due to the lithium ions coordinated to FSI anion, but occurring at different radial distance from the center of mass of the anion since the coordination sites are not equivalent with respect to the center of mass. Focusing on the CPME-based electrolytes where the peaks below 500 pm are more defined, the first peak at ∼290 pm corresponds to Li⁺ coordinated to FSI⁻ O atoms bonded to different S atoms (Fig. S1C). The second peak at ∼390 pm is due Li⁺ coordinated to only one oxygen atom (Fig. S1D). The RDFs for DME based electrolytes on the other hand tell a different story. For solutions with less than about 3.5 M LiFSI, the RDFs lack any peaks below 500 pm. At 3.5 M, the peaks gradually increase as the concentration of LiFSI increases, and at the same time the RDF peaks at longer distances decrease. This clearly shows that different modes of interaction exist at low and higher concentrations of LiFSI in DME. It is noteworthy that simulations with standard non-scaled charges for ions, the peaks below 500 pm persist at all concentrations of LiFSI in DME considered in this study, which is against experimental reality (see SI Fig. S2). This is because interactions between solvent and ions are overestimated in non-scaled charges which don't take into account solvent polarization. The total absence of the peak in scaled charges at low concentrations indicates that the FSI anion and Li⁺ exist as ‘free’ ions solvated by DME in this case. The center of mass RDFs between FSI⁻ and the solvents at different concentrations are also shown in Fig. 1. These also show significant differences in the case of DME electrolytes at different concentrations. The Li/solvent RDFs do not show any significant variations at different concentrations (Fig. S3). The coordination numbers of different species around the cation and ion are shown in Table 1 for both types of electrolytes at 1 M, 2.5 and 5 M concentrations. These were calculated from based on the first minima of the radial distribution functions. The calculated solvent coordination numbers around lithium ion are in good agreement with those calculated from Raman spectra.^9,37 The coordination numbers of Li around FSI anion also qualitatively agree with experimental data, in that CPME based electrolyte forms aggregates and counterion pairs at lower concentrations than DME based electrolyte.⁸

Fig. 1 Radial distribution functions of different pairs at different concentrations: FSI⁻ (centre of mass)/Li⁺ in CPME (A), FSI⁻ (centre of mass)/Li⁺ in DME (B), FSI⁻ (centre of mass)/CPME (centre of mass) (C) and FSI⁻ (centre of mass)/DME (centre of mass) (D).

Table 1 Calculated coordination numbers for different pairs in LiFSI/CPME and LiFSI/DME electrolytes at different concentrations

	Li^+-solvent			FSI^−-solvent			FSI^−-Li^+
	100/900	250/750	500/500	100/900	250/750	500/500	100/900	250/750	500/500
CPME	2.44	1.94	0.98	12.65	10.74	6.86	0.79	1.53	2.84
DME	2.10	2.05	1.00	17.50	18.83	8.53	—	—	2.2

---

FSI⁻/Li⁺ aggregation analysis

To further compare with experimental data, cluster size distribution was calculated using Travis software.^38,39 In this case the reference molecule is FSI⁻ anion. To do this, 5000 reference molecules are selected randomly from the available trajectory frames, ensuring a wide sampling strategy. The selection is done before the analysis to remove bias. The solvation sheath of each reference molecule is then analyzed based on the first coordination shells of FSI⁻/Li⁺ and FSI⁻/solvent radial distributions functions. The resulting distribution of clusters of different forms of FSI⁻ are shown in Fig. 2. Note that CIP refers to counter ion pairs, AGG-I refers to FSI⁻ coordinated to 2 Li⁺ ions and AGG-II refers to FSI⁻ surrounded by more than 2 Li⁺ ions. The one Li⁺ per FSI⁻ O atom (Fig. S1D) is more favorable for aggregates formation, hence the observed increase in the second RDF peak as concentration increases. The DME-based electrolyte is dominated by Free FSI⁻ at low concentrations in agreement with experimental data. In the case of CPME-based electrolytes, CIPs are dominant at low concentrations with significant amounts of AGG-I and AGG-II. Experimentally, aggregates are dominant even at low electrolyte concentrations for CPME. Given a simplified treatment of polarization in the simulation, we consider this satisfactory agreement with experimental data.

Fig. 2 FSI⁻ species distribution in different concentrations of LiFSI in CPME (A) and DME (B). FSI⁻ can appear as free and not bound to Li⁺, bound to 1 Li⁺ (CIP), bound to 2 Li⁺ (AGG-I) or bound to 3 or Li⁺ ions (AGG-II).

Solvent distribution and aggregation

The change in the FSI⁻/DME center of mass RDFs from low to high concentration suggests that the interaction between FSI⁻ and DME could have an impact in the formation of FSI⁻/Li⁺ aggregates. It is reasonable to expect that in the case of strong FSI/solvent interactions, the solvent can compete with Li⁺ for FSI⁻ interaction sites. Another possibility is that FSI⁻/Li⁺ aggregates may be stabilized by solvents molecules via secondary interactions of the FSI⁻ coordinated Li⁺ with a solvent molecule. In order to investigate both possibilities, the distributions of solvent molecules around the CIPs and aggregates were calculated as above from MD trajectories and these are shown in Fig. 3. Note the relatively small number of CIPs formed in 1 M LiFSI in DME compared to other systems, again indicating that most of the ions exist as free ions. The CIPs in 1 M LiFSI in DME are mostly solvated by 2 DME molecules. For 1 M LiFSI in CPME and 5 M LiFSI in DME, the CIPs and aggregates are mostly solvated by 3 solvent molecules. The 5 M LiFSI in CPME solution has both CIPs and aggregates solvated most likely by 2 CPME molecules. It is difficult to conclude from these results if the solvent plays any role in stabilizing the aggregates. The fact that for one concentration and solvent type, the most probable number of solvents around CIP and aggregates is the same suggests minimal solvent role once these are formed. If solvent stabilization played a significant role, one would expect that aggregates would be bound to more solvent molecules than CIPs.

Fig. 3 Distribution of solvent molecules around (Li⁺)_nFSI⁻ counter ion pairs (n = 1), aggregates (n > 1) and free FSI⁻ (n = 0) in the first solvation shell in 0.1 M LiFSI in DME (A), 0.1 M LiFSI in CPME (B), 0.5 M LiFSI in DME (C) and 0.5 M LiFSI in CPME (D).

There are more DME molecules around FSI⁻ possibly because DME is less sterically hindered and more flexible than CPME. The better flexibility of DME is evidenced by the dihedral distribution functions of the C(H₃)–O–C–C and C–O–O–C in CPME and DME shown in Fig. 4. For DME, significant differences are noted in the dihedral angle at different concentrations. Unlike DME, the distribution for CPME show very little changes with concentration.

Fig. 4 Dihedral angle distribution functions: C(H₃)–O–C–C in CPME (A) and C–O–O–C in DME (B).

FSI⁻/Li⁺ and Li⁺/solvent interactions

Since the interaction between solvent and anion are not significant to lead to the observed differences between CPME and DME electrolytes interactions, binding energies (BE) for FSI⁻/Li⁺ and Li⁺/solvent species were calculated from the energies of quantum mechanically optimized geometries using the equation

BE = E_complex − ∑E_monomers (2)

Here E_complex refers to the energy on the optimized bound entities and E_monomers are the energies of the individually optimized units.⁶ Since the interactions between the species are expected to be weaker in solution, implicit solvation was included via the Conductor-like Polarizable Continuum Model (CPCM).⁴⁰ The calculated binding energies are shown in Fig. 5.

Fig. 5 Calculated binding energies in gas phase and in solution using CPCM implicit solvation model: (A) Li⁺-(CPME)_1–3 and Li⁺-FSI⁻ in CPME, (B) Li⁺-(DME)_1–3 and Li⁺-FSI⁻ in DME.

The binding energies for Li⁺-DME and Li⁺-CPME are in close agreement with previous calculations.^9,37 The binding between FSI⁻ and Li⁺ in gas phase is quite strong, −137.5 kcal mol⁻¹, but is severely weakened in solution. The binding is stronger in CPME than in DME. A more interesting observation though is that the binding between FSI⁻ and Li⁺ in DME solution is almost equal to that of Li⁺-DME in solution and much weaker than in Li⁺-(DME)₂ and Li⁺-(DME)₃. Li⁺-(CPME)₂ and Li⁺-(CPME)₃ also bind strongly than FSI⁻/Li⁺ in solution. However, analysis of the DME/Li⁺ ion solvation in molecular dynamics trajectories at 1.0 M shows that the clustering is predominantly Li⁺-(DME)₂ (90%), the remainder being Li⁺-(DME)₃. This is consistent with the average coordination number of 2.1 shown in Table 1. In contrast, for the CPME based electrolyte solution at the same concentration, Li⁺-(CPME)₂ and Li⁺-(CPME)₃ account for 54% and 45% respectively (see SI Fig. S4). Non-solvated Li⁺ only accounts for 1%. Again, this is consistent with the coordination number of 2.44 shown in Table 1. Previous simulations of 1.0 M LiPF₆ in DME showed that the Li⁺ solvation is completely composed of Li⁺-(DME)₃.⁶ The study did not consider polarization effects as done in the current study. Given the different electrostatics between Li⁺/FSI⁻ and Li⁺/PF₆⁻, differences in clusters size distribution is still expected. At higher concentrations more aggregates are formed for both CPME and DME based electrolytes, however, the solvation shells in DME based electrolytes still contain more solvent molecules than in CPME. For battery the CPME based electrolyte is still more desirable as involvement of less solvent molecules in the solvation shell can lead to a more stable SEI.⁴¹

Statistically, higher number of solvent molecules around a Li⁺ ion is beneficial to the formation of aggregates especially for dilute solutions. This is because the “lithium demand” by solvent is less for electrolytes that are dominated by higher number of solvent molecules around lithium solvation shell. Low demand for Li⁺ by solvent molecules means that more Li⁺ ions are available to form aggregates and counter ion pairs with FSI⁻ anion. Beside the fact that binding between Li⁺ and FSI⁻ in DME is weaker than binding in Li⁺-(DME)₂, the predominance of Li⁺-(DME)₂ could explain why the DME based electrolyte lacks the first Li solvation shell at low concentrations (Fig. 1). The Li⁺-(DME)₂ clusters with four coordinate Li⁺ ion are so stable that more Li⁺ ions must be available to satisfy the number of available DME molecules. In this case the ratio LiFSI : DME must be greater than 1 : 2. From the radial distribution functions in Fig. 1, this condition is only satisfied for concentrations of LiFSI in DME above 0.35 M. Indeed, it is only at concentrations above this that the first solvation shell begins to appear.

It is possible to gauge the stability of the clusters in solutions by tracking the Li⁺ ions and solvent molecule. Indexing each Li⁺ and solvent molecule makes it possible to know how long a selected Li⁺ ion maintains a unique solvation environment. Fig. 6 shows the average survival times of all the unique Li⁺-solvent, Li⁺-(solvent)₂, Li⁺-(solvent)₃ and Li⁺-FSI⁻ in 1.0 M LiFSI electrolytes. As expected, Li⁺-(DME)₂ clusters survive much longer than Li⁺-(DME)₃ clusters. In the case of CPME, the Li⁺-(CPME)₂ clusters survive slightly longer than Li⁺-(CPME)₃ clusters. These observations also explain why the aggregation still occurs in CPME based electrolytes at low concentration even though binding between Li⁺ and FSI⁻ is weaker than in Li⁺-(CPME)₃; the latter has a much shorter life compared to Li⁺-FSI⁻. Since the solutions have different viscosities, comparison is made only between clusters within the same electrolyte.

Fig. 6 Average survival times in picoseconds for solvated Li⁺ clusters and Li⁺-FSI⁻ in different DME and CPME.

Conclusions

The interactions between Li⁺ ion, bis(flurosulfonyl)imide anion (FSI⁻) and solvent molecules in cyclopentylmethyl ether (CPME) and 1,2-dimethoxy ethane (DME) have been studied using molecular dynamics and ab initio simulations. The results show that inclusion of the polarization effects in force field based molecular dynamics and implicit solvation in ab initio simulations leads to a good agreement with experimental data. Clustering analysis shows that at low concentrations of LiFSI in DME, the Li⁺ predominantly exists as a free ion solvated by two DME molecules (90% Li⁺-(DME)₂) and the binding in this system is not only much stronger than between Li⁺ and FSI⁻ in the same solution, but the cluster itself survives for a longer time. In contrast, in CPME based electrolytes both Li⁺-(CPME)₂ and Li⁺-(CPME)₃ are present in almost the same amounts, but in this case the binding energy of Li⁺-FSI⁻ is about the same as Li⁺-(CPME)₂, while Li⁺-(CPME)₃ has a very short survival lifetime. The onset concentration of the formation of aggregates in DME based electrolytes suggests that all the DME molecules existing in solution in principle should be coordinated to Li⁺ and only the remaining Li⁺ are then free to bind to FSI⁻ to form aggregates and counter ion pairs.

Author contributions

Conceptualization, F. M.; methodology, F. M; validation, F. M and M. N.; formal analysis, F. M. and M. N.; investigation, F. M. and M. N.; writing – original draft preparation, F. M. and M. N.; writing – review and editing, visualization, F. M. and M. N. Both authors have read and agreed to the published version of the manuscript.

Conflicts of interest

We wish to confirm that there are no known conflicts of interest associated with this publication.

Data availability

The datasets generated and analyzed during the current study are available within this published article and its supplementary information (SI). Supplementary information: these include geometrical structures of the investigated solvents, selected radial distribution functions of simulations without scaled charges and Li⁺-(SOLVENT)n (n = 1–4) distribution in both DME and CPME based electrolytes. See DOI: https://doi.org/10.1039/d5ra08464d.

Acknowledgements

The authors acknowledge the support of Botswana International University of Science and Technology.

References

1. R. Younesi, G. M. Veith, P. Johansson, K. Edström and T. Vegge, Energy Environ. Sci., 2015, 8, 1905–1922.
2. X. Tang, S. Lv, K. Jiang, G. Zhou and X. Liu, J. Power Sources, 2022, 542, 231792.
3. Q. Li, Z. Cao, W. Wahyudi, G. Liu, G.-T. Park, L. Cavallo, T. D. Anthopoulos, L. Wang, Y.-K. Sun, H. N. Alshareef and J. Ming, ACS Energy Lett., 2021, 6, 69–78.
4. U. S. Meda, L. Lal, S. M and P. Garg, J. Energy Storage, 2022, 47, 103564.
5. G. Li, X. Lyu, A. Nguyen, R. Kou, C. George, S. Wu, R. Li, K. Wang, T. Li and D. Wang, Adv. Energy Mater., 2025, 15(19), 2405680.
6. X. Chen, N. Yao, B. S. Zeng and Q. Zhang, Fundam. Res., 2021, 1, 393–398.
7. Z. Wang and B. Zhang, Energy Mater. Devices, 2023, 1(1), 9370003.
8. H. V. Ramasamy, S. Kim, E. J. Adams, H. Rao and V. G. Pol, Chem. Commun., 2022, 58, 5124–5127.
9. H. Zhang, Z. Zeng, F. Ma, Q. Wu, X. Wang, S. Cheng and J. Xie, Angew. Chem., 2023, 62(21), e202300771.
10. Z. Wang, S. Weng, H. Zhang, L. Wang, X. Yao, H. Tu, D. Huang, S. Lu, L. Liu, J. Xue, F. Zhang, G. Wu, J. Zheng, Q. Wang, L. Chen, J. Xu, H. Li and X. Wu, Adv. Sci., 2025, 12(37), e17305.
11. K. Watanabe, N. Yamagiwa and Y. Torisawa, Org. Process Res. Dev., 2007, 11, 251–258.
12. E. Park, J. Park, K. Lee, Y. Zhao, T. Zhou, G. Park, M.-G. Jeong, M. Choi, D.-J. Yoo, H.-G. Jung, A. Coskun and J. W. Choi, ACS Energy Lett., 2023, 8, 179–188.
13. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
14. S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. Van Der Spoel, B. Hess and E. Lindahl, Bioinformatics, 2013, 29, 845–854.
15. L. S. Dodda, I. Cabeza de Vaca, J. Tirado-Rives and W. L. Jorgensen, Nucleic Acids Res., 2017, 45, W331–W336.
16. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
17. L. S. Dodda, J. Z. Vilseck, J. Tirado-Rives and W. L. Jorgensen, J. Phys. Chem. B, 2017, 121, 3864–3870.
18. J. N. Canongia Lopes and A. A. H. Pádua, J. Phys. Chem. B, 2004, 108, 16893–16898.
19. A. S. L. Gouveia, C. E. S. Bernardes, L. C. Tomé, E. I. Lozinskaya, Y. S. Vygodskii, A. S. Shaplov, J. N. C. Lopes and I. M. Marrucho, Phys. Chem. Chem. Phys., 2017, 19, 29617–29624.
20. J. N. Canongia Lopes and A. A. H. Pádua, Theor. Chem. Acc., 2012, 131, 1129.
21. L. X. Dang, J. Chem. Phys., 1992, 96, 6970–6977.
22. Z. R. Kann and J. L. Skinner, J. Chem. Phys., 2014, 141(10), 104507.
23. I. V. Leontyev and A. A. Stuchebrukhov, J. Chem. Theory Comput., 2010, 6, 3153–3161.
24. I. V. Leontyev and A. A. Stuchebrukhov, J. Chem. Theory Comput., 2012, 8, 3207–3216.
25. C. Park, M. Kanduč, R. Chudoba, A. Ronneburg, S. Risse, M. Ballauff and J. Dzubiella, J. Power Sources, 2018, 373, 70–78.
26. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
27. H. J. C. Berendsen, in Computer Simulation in Materials Science, Springer Netherlands, Dordrecht, 1991, pp. 139–155.
28. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126(1), 014101.
29. M. A. Cuendet and W. F. van Gunsteren, J. Chem. Phys., 2007, 127(18), 184102.
30. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
31. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
32. F. Neese, J. Comput. Chem., 2023, 44, 381–396.
33. B. Helmich-Paris, B. de Souza, F. Neese and R. Izsák, J. Chem. Phys., 2021, 155(10), 104109.
34. F. Neese, F. Wennmohs, A. Hansen and U. Becker, Chem. Phys., 2009, 356, 98–109.
35. F. Neese, J. Comput. Chem., 2003, 24, 1740–1747.
36. Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 120, 215–241.
37. G. Zhang, T. Zhang, Z. Zhang, R. He, Q. Wang, S.-S. Chi, Y. Cui, M. D. Gu, Z. Liu, J. Chang, C. Wang, K. Xu and Y. Deng, Nat. Commun., 2025, 16, 4722.
38. M. Brehm, M. Thomas, S. Gehrke and B. Kirchner, J. Chem. Phys., 2020, 152(16), 164105.
39. M. Brehm and B. Kirchner, J. Chem. Inf. Model., 2011, 51, 2007–2023.
40. V. Barone and M. Cossi, J. Phys. Chem. A, 1998, 102, 1995–2001.
41. Q. Wang, J. Wang, J. R. Heringa, X. Bai and M. Wagemaker, ACS Energy Lett., 2024, 9, 3796–3806.

---