Improved modeling of battery electrolytes: betting on model fitting or quantum effects?

Abstract

The accurate modeling of solvent dynamics and ionic interactions is of crucial importance for the development of novel electrolytes in next-generation metal-ion batteries. This study presents a critical evaluation of the semi-classical computational approach, the adaptive quantum thermal bath (adQTB) method, as a methodology for capturing the key properties of glyme-based solvents and their Ca²⁺-based electrolyte solutions. Simulations reveal that the adQTB method is particularly effective in accurately reproducing vibrational spectra, while offering good transferability across systems and conditions without requiring empirical parameter adjustments. In the context of electrolyte solutions, semi-classical adQTB simulations in combination with graph theory analysis indicate a distribution of the various charge-carrying clusters that is closely aligned with the conductivity measurements previously reported [Nguyen et al., Phys. Chem. Chem. Phys., 2022, 24, 21601], in sharp contrast to the empirically scaled force field. These findings emphasize the necessity of incorporating nuclear quantum effects for reliable electrolyte modeling, thereby paving the way for the advancement of post-Li battery technology.

---

In the field of metal and metal-ion batteries, the development of liquid, polymer, or solid electrolytes that address both safety and lifespan concerns has emerged as a significant area of research.^1–5 Over the past decade, this field has experienced substantial growth, driven primarily by advancements in post-Li technologies that offer promising alternatives to conventional Li-ion batteries.^6–9 As the demand for more efficient and safer energy storage solutions continues to increase, a comprehensive understanding of the mechanisms of ionic conductivity and the influence of solvation properties on the possible degradation of electrolytes during electrochemical charge/discharge cycles has become critically important. Recent studies underscore the need to model electrolyte behavior under diverse conditions, including varying concentrations, temperatures, and external potentials.^10–13 While ab initio molecular dynamics (AIMD) is widely regarded as the gold standard, its computational cost makes it untractable for such high-throughput screening. Classical molecular dynamics (MD) then becomes the method of choice, provided that the simulation setup used to investigate the interactions and the dynamic behavior of ions and solvent molecules at the atomic scale accurately capture the relevant physical effects. These simulations offer atomistic insights into key phenomena such as ion transport,^14–17 solvation dynamics,^18–20 and interfacial processes,^21–23 all of which are critical to battery performance.^24–30 Classical MD simulations thus serve as a virtual laboratory, offering a cost-effective means to complement experimental studies and improve our understanding of metal-ion batteries.

In this context, the development of accurate force fields to capture the dynamics of electrolyte solutions remains a vital and ongoing area of research.^31–39 Polarizable force fields such as Drude oscillator,^40–45 or induced dipole models^46–50 are required to explicitly account for the atomic polarizability and to capture the molecular response to various electric fields, including static and dynamic effects. This is particularly relevant for electrolytes based on low-dielectric organic solvents, where long-range electrostatic interactions dominate ion transport and solvation phenomena.¹¹ In a recent work, we investigated the glyme solvents (glycol methyl ethers), a class of low dielectric solvents used in post-Li battery electrolytes.^51,52 Using the AMOEBA model, we introduced a transferable fitting procedure to develop a universal interaction potential for the entire glyme family (–O–(–CH₂–CH₂–O–)_n–) from small oligomers (Gn) to the PEO polymer (G∞ = PolyEthylene Oxide). Results showed that the AMOEBA force field well captures local bonding interactions and solvation structures, while unexpectedly overestimating macroscopic properties, such as density and dielectric constant. This discrepancy underscored the need for further refinement to improve the alignment of computational predictions with experimental data and our understanding of electrolyte behavior at both molecular and macroscopic scales.

Accurately capturing nuclear quantum effects (NQE), both in ab initio or classical MD simulations is a critical step when modeling light atoms, like hydrogen, especially at low temperatures.^53–57 Among existing methodologies that explicitly incorporate NQE, path integral molecular dynamics (PIMD) stands out as a powerful, albeit computationally expensive, approach that accounts for all possible quantum paths of the system. A widely used approximation of PIMD is ring polymer molecular dynamics (RPMD),⁵⁸ which makes use of multiple interacting copies of the system. However, the computational cost of RPMD still remains prohibitive for large systems, particularly at low temperatures where the number of required copies increases with the temperature decrease. This issue has driven the search for more practical alternatives. The quantum thermal bath (QTB) method,⁵⁹ based on Langevin dynamics, offers a computationally affordable solution with limited extra-cost compared to NQE-free MD simulations. Unlike the classical Langevin thermostat, which applies random forces (white noise) on all nuclei,⁶⁰ the QTB method incorporates frequency-specific adjustments via colored noise to ensure a reliable (i.e. quantum-mechanics) energy partition of vibrational modes. By modifying the random-force amplitudes, the QTB method can supply or withdraw quanta of energy from each mode, aiming to reproduce the correct zero-point-energy (ZPE). However, the original QTB approach, which treats vibrational modes in a unified manner, leads to the so-called ZPE leakage, particularly from high to low-frequency modes,⁶¹ thereby compromising the method's accuracy and applicability. To address this challenge, Mangaud et al.⁶² introduced the adaptive quantum thermal bath (adQTB), an advanced version of the QTB method. The adQTB incorporates frequency-dependent friction coefficients, derived from the fluctuation-dissipation theorem, which are dynamically adjusted (i.e. “on the fly”) during the MD simulations to prevent ZPE leakage. By overcoming this limitation, adQTB significantly improves upon traditional QTB, therefore offering better accuracy. As a result, it represents a valuable tool for incorporating nuclear quantum effects into MD simulation of complex electrolyte systems, with no substantial extra-cost.

In this letter, we discuss the influence of explicitly incorporating NQE via the adaptive quantum thermal bath (adQTB) for two systems having experimental reference data available: pure liquid glymes^63–66 and glyme-based electrolyte solutions with a Ca(TFSI)₂ salt.¹¹ Results are compared to the commonly used van der Waals (vdW) scaling approach,^67–69 which adjusts the vdW parameters of the force field to fit various experimental data. Although the vdW-scaling approach is popular and relatively straightforward to implement, the parameter optimization process relies on the availability and accuracy of experimental data. Moreover, such empirically-optimized force fields often lack transferability across different systems, hence limiting their versatility compared to methods introducing the quantization of nuclei either from first-principles (e.g. PIMD) or in an effective manner (e.g. adQTB). Despite its limitations, the vdW-scaling approach can still capture key electrolyte properties, as nuclear quantum effects are inherently reflected in the experimental data used to fit the force field parameters. While both vdW-scaling and adQTB approaches are expected to improve the match between simulations and experiments compared to unscaled and purely classical approaches, the correction offered by each method may vary, particularly at low temperatures where nuclear quantum effects are significant. The Boltzmann distribution on its own becomes a more satisfactory description of the energy partition at higher temperatures, where the gap between classical and semi-classical approaches should close.

For the sake of clarity, the previously reported AMOEBA force field will be used as a reference and hereafter referred to as “initial”. The force field resulting from the scaling of its vdW parameters will be denoted “vdW-sc”, while the use of adQTB with the initial force field will be denoted “adQTB”. The performance of these three approaches is exemplified in the following on local and macroscopic properties, including bond distributions, vibrational structure, density, dielectric constant, and solvation structures. For this specific study, a vdW-scaling factor of k = 1.06 was extracted from the experimental tetraglyme (G4) densities,⁶⁴ as the estimation of this quantity is highly sensitive to intermolecular interactions (see SI, Fig. S1). As shown in Fig. 1, this scaling factor improves the reproduction of the experimental density for both diglyme (G2) and tetraglyme (G4) compared to the initial approach, with relative errors below 1.4% across the temperature range from 273 to 363 K. Incorporating nuclear quantum effects also improves the results, though densities are still slightly overestimated (around 2% error for both G2 and G4). As expected, the adQTB method offers a good transferability across systems and temperatures with relative errors remaining stable in the whole test-case study (see SI, Tables S1 and S2 for exact density values).

Fig. 1 Evolution of the liquid diglyme G2 (top) and tetraglyme G4 (bottom) densities with respect to temperature, as obtained from MD simulations using the initial (yellow), vdW-sc (green) and adQTB approaches (red). The horizontal dotted lines represent the experimental density⁶⁴ for both systems at different temperatures. Relative errors compared to experiments are also indicated on each plot.

To have a more complete comparison of the three approaches, it is interesting to study the macroscopic properties that were not explicitly considered for the vdW-parameter optimization. As shown in Table 1, both vdW-sc and adQTB perform similarly for the vaporization enthalpies of G2 and G4, and slightly improve the initial results. For the dielectric constants, results are consistent with the density results, i.e. the larger the density, the larger the dielectric constant. This is consistent with the positive linear relationship between density and dielectric constant (illustrated in Fig. S2) as both properties vary with inverse volume. Compared to the density results, the larger differences for the dielectric constants observed between the different methods suggest that dipole moment fluctuations are highly sensitive to even slight changes in intermolecular interactions. Note that the good agreement between experimental and simulation data for the initial valence-fit force field without accounting for NQE (both for G2 and G4) reconcile the conflicting results discussed in our previous work⁵¹ where dielectric constants were erroneously overestimated.⁵²

Table 1 Vaporization enthalpies (ΔH_vap in kJ mol⁻¹) and dielectric constants (ε) calculated at 298.15 K for liquid G2 and G4 using the initial, vdW-sc and adQTB approaches and compared to experimental values^65,66

	Exp.	Initial	vdW-sc	adQTB
ΔHvap,G2	54.6 ± 3.8%	47.3 ± 8.0%	51.4 ± 8.2%	54.4 ± 9.6%
ΔHvap,G4	76.9 ± 3.4%	90.2 ± 2.6%	89.3 ± 2.9%	87.6 ± 1.8%
εG2	7.40 ± 0.10	7.34 ± 0.03	6.10 ± 0.19	6.76 ± 0.03
εG4	7.78 ± 0.01	8.10 ± 0.25	6.34 ± 0.15	7.98 ± 1.15

---

Quantifying nuclear motion is crucial for vibrational properties, especially for intramolecular vibrations. The difference between the classical and semi-classical methods is reflected in the simulated infrared (IR) spectra of liquid G2 and G4 given on Fig. 2. The discrepancy between computed and experimental relative intensities reflects an intrinsic limitation of the force field approach, which cannot reproduce charge redistribution along vibrational modes. This may result in inaccurate dipole moments variations as illustrated in SI, Tables S4 and S5 for the C–H bonds. The initial and vdW-sc approaches lead to nearly superimposed IR spectra in the low- and high-frequency regions, consistent with the fact that the vdW-scaling primarily alters intermolecular interactions (see SI, Fig. S3, S4 and Table S3). The adQTB method shows clear improvement in the prediction of peaks frequency and broadening with a red-shift of the high-frequency modes around 1110, 1455, 2880 cm⁻¹ for both G2 and G4. This red-shift is particularly significant for the C–H stretching modes around 2880 cm⁻¹, which lower by approximately 100 cm⁻¹, hence falling within the expected range of 2750–3000 cm⁻¹.

Fig. 2 IR spectra calculated for the liquid diglyme (G2) (top) and tetraglyme (G4) (bottom) as explained in computational details using the vdW-sc (green) and adQTB (red) methods. The dotted lines correspond to experimental peaks⁶³ (see SI, Fig. S3 and S4 for full experimental spectra, and Table S3 for peak lists).

Overall, these results underline the importance of a semi-classical approach for accurately capturing the vibrational behavior of light atoms such as hydrogen. This improved vibrational representation is associated with the adQTB's ability to sample a wider range of bond lengths, as demonstrated in Fig. 3 for both G2 and G4 chain lengths. Such a wider sampling suggests that incorporating NQE should be particularly beneficial for highly hydrogenated organic compounds used as electrolyte solvents in post-Li batteries or any other organic systems, including ether chains, hydroxyl groups, alkyl chains, and even cyclic alkanes.

Fig. 3 Distribution of intramolecular C–H distances for all CH₂ and CH₃ groups of diglyme (left) and tetraglyme (right) molecules, using the vdW-sc (green) and adQTB (red) methods.

Another interesting question is whether the changes observed in the average response of molecular bonds due to the method specificity have any impact on the overall conformation of individual glyme chains. In the case of the highly flexible tetraglyme (G4), a key intramolecular parameter is the extent to which its conformation tends toward linearity or crown shape. As shown in Fig. 4, the overall profile of the density graph remains almost identical for all three simulation setups with, however, a slight elongation of the most abundant CH₃⋯CH₃ distances for vdW-sc. This corresponds to a small linearization of the ether chains which can be more likely correlated to the slightly lower density obtained with this force field (see Fig. 1), than to a fundamental difference in the conformational exploration during the MD simulation.

Fig. 4 Distribution of intramolecular end group (CH₃⋯CH₃) distances extracted from room-temperature MD simulations of bulk liquid using the initial (yellow), vdW-sc (green) and adQTB (red) approaches.

In the context of electrochemical energy storage, the investigation of the solvation environment of charge carriers is particularly compelling, as a change in dominant structures may affect ion diffusion within battery electrolytes.^11,70 The nature and the relative occurrence of these charge carriers may depend on the force field and on whether NQE are explicitly taken into account in the MD simulations. For this purpose, an electrolyte model of ∼1 M of Ca(TFSI)₂ salt in tetraglyme (G4) was studied at different temperatures, with the initial, vdW-sc and adQTB methods (see SI, Section S3 for TFSI⁻ force field parameters, and Table S7 for exact molar concentrations). A graph theory based code, developed by Vatin et al.,⁷¹ was then used to identify and count the different contact and solvent-separated ion pairs occurring in the solution, in order to determine the distribution of neutral and charged species. Based on the radial distribution functions (RDF), different cutoff distances characteristic of the cation–anion (Ca²⁺⋯O(TFSI⁻)), cation–solvent (Ca²⁺⋯O(G4)) and solvent–anion (O(G4)⋯O(TFSI⁻)) interactions were applied to capture the different pair associations that may occur in the solution (see SI, Fig. S5 and Table S6). The charge carrying clusters containing at least one Ca²⁺ ion were then reconstructed, and their total charges were deduced from the presence of one or more TFSI⁻ anion as contact or solvent-separated counter-ions.

As illustrated in Fig. 5, the relative abundance of the charged and neutral clusters strongly depends on the force field and the incorporation of NQE. Indeed, noticeable differences are observed between the different approaches, not only in the species relative ratio, but also in their evolution upon temperature increase. The vdW-sc approach clearly promotes ion-pair association, as evidenced by the constant decrease of 2+ charged clusters as temperature increases. The distribution obtained with the initial method qualitatively compares with vdW-sc with, however, a smoother decrease in the number of 2+ charged clusters with temperature increase. For both methods, this monotonic temperature-dependent ion-pairing behavior leads to a continuous increase in the number of neutral clusters, with neutral species becoming dominant at the highest simulated temperature. This is in sharp contrast with the adQTB approach which first displays a notable augmentation of 2+ charged species with temperature increase followed by a diminution above a temperature threshold between 328 and 388 K. At 298 K, the average cluster size and the occurrence of multi-cationic clusters increase greatly by the inclusion of NQEs, highlighting their influence on long-range interactions (see Table S8 and Fig. S6). Between 298 K to 328 K, the adQTB also shows a 14% increase in free TFSI⁻ anions in the electrolyte, suggesting a more complex dependence of cluster formation on temperature than the mere acceleration of ion pairing kinetics. This is consistent with the non-monotonous variation of the average cluster size and the occurrence of multi-cationic species (Table S8) with temperature across all three approaches.

Fig. 5 Distribution of charges carried by the clusters reconstructed by graph theory from MD simulations performed at different temperatures. From top to bottom, initial force field (yellow), vdW-scaled force field (green), initial force field with adQTB (red).

According to the experimental data reported by Nguyen et al.,¹¹ the ionic conductivity of Ca(TFSI)₂@G4 electrolyte solutions first increases with temperature, prior to slightly decreasing above a temperature threshold which rises with the salt concentration. For concentrations ranging from 0.5 to 0.8 M, this temperature threshold is consistently observed between 343, 383 K which is consistent with the adQTB results. The amount of charge carriers being related to the ionic conductivity, all other parameters assumed equal for a given temperature, the trend obtained with the vdW-sc approach – and to a lower extent with the initial method – deviates from experimental data, therefore indicating a substantial interdependence between the vdW and other intermolecular parameters of the AMOEBA model. This indicates that the application of a scaling factor to selected-only contributions has a detrimental impact on the accuracy and transferability of the force field and cannot fully replace the explicit inclusion of nuclear quantum effects in complex systems like electrolytes. This is consistent with recent reports showing that the explicit inclusion of NQE, using Path Integral MD method, improves self-diffusion constant estimation.⁷² Although the current implementation of adQTB is not suitable for the calculation of dynamic properties such as conductivity, these results suggests that indirect ionic conductivity predictions may be practical for various electrolyte solutions, at least from a comparison perspective, based on the statistics of charged species distributions. It is known that different ion association motifs, ranging from contact or solvent separated ion pairs to extended clusters, behave differently in the electrolyte, affecting the transport properties as well as the possible decomposition mechanisms.^73–76 As briefly illustrated in Fig. S7, we have a zoology of clusters that carry the same number of charge. Their distinction and identification via graph theory based methods^71,77,78 is therefore of interest for future work. It is not possible to directly associate these clusters to diffusion constants as suggested by previous works,^73,74 due to the limitations of the adQTB framework. However, one could aim to approximate diffusion constants and ionic conductivity by exploiting the relevant physical properties of the clusters such as mass, charge, and morphology, as well as those associated to their environment such as viscosity and temperature. Machine-learning algorithms may prove invaluable in this respect to predict trends in conductivity as a function of these properties. Such an approach could provide an efficient, albeit qualitative, way to screen novel electrolytes currently under development for Li and post-Li batteries, without engaging in computationally expensive simulations of transport properties.¹⁴

In conclusion, this study compares the direct incorporation of nuclear quantum effects (NQE) via the adaptive quantum thermal bath (adQTB) method with the conventional approach of scaling van der Waals (vdW) parameters, evaluating their influence across three key areas: macroscopic properties, spectroscopic response, and molecular conformations. For pure liquids, both methods performed comparably on a macroscopic level. However, adQTB offers a critical advantage in transferability, as it does not require system-specific parameter adjustments, whereas vdW-scaling does. This is especially valuable for systems lacking experimental data, such as those at an early stage of development. Nonetheless, for chemically similar systems, such as different glymes, a predetermined scaling factor remains valid for a specific macroscopic property. On a molecular level, vdW-scaling alone proves inadequate for accurately capturing high-frequency vibrational peaks in the infrared spectrum, where the explicit inclusion of NQE significantly improves peak frequency and broadening with respect to experimental spectra.

For chemically diverse systems such as electrolyte solutions, the vdW-scaling approach may fall short, as the vdW-scaled force field predicts significantly different solvation environments for charged species at different temperatures in a way that contradicts experimental observations. In contrast, the adQTB method is in line with experimental trends and could be used as a powerful way to qualitatively estimate ionic conductivity with respect to increasing temperature while predicting different occurrences of charged species. The influence of NQE accounted for through adQTB therefore goes beyond a simple shift in simulated temperature that could be assimilated to the zero point energy. Further investigations are needed to understand how NQE impact long-range interactions in the electrolyte, particularly with respect to charge distribution among supermolecular carriers. In particular, the accurate determination of the distribution of various aggregates and charges could be a particularly fruitful endeavour. Coupled with machine-learning algorithms, the characteristic descriptors of these aggregates (size, charge, composition) determined by graph theory could significantly enrich the databases, with the aim of moving towards more rational design of electrolytic solutions.

Computational details

Molecular dynamics simulations were performed using the Tinker-HP program package (version 1.2).⁷⁹ Liquid densities were calculated by performing 1 ns of NPT simulation and averaging the value over the last 125 ps of the trajectory. It has been verified that this average remains within 2 standard deviations (σ) of the previous 250 ps. A time step of 1 fs has been used.

An NVT simulation of 48 ns was then performed on the NPT equilibrated system, conserving the 1 fs time step. The dielectric constants were calculated over the two 24 ns halves, based on the fluctuation of the total dipole moment as shown in eqn (1):

where for polarizable models such as AMOEBA, the optical dielectric constant ε_∞ is given by the Clausius–Mossotti relation:

The molecular polarisability α_i was calculated using the polarize module of Tinker (version 8.10).⁸⁰ In the case of G4, the same NVT trajectory was used to calculate the average potential energy of liquid (U_liq) for the vaporization enthalpy (ΔH_vap,G4). A 4 ns NVT trajectory of a single G4 molecule in a 50 × 50 × 50 Å box was performed to obtain the potential energy for the ideal gas (U_gas). The average values for U_gas and U_liq were calculated over the last 1/4 of their respective trajectories, ensuring that they stay within 2σ of the previous average.

As suggested by Wang & Hou⁸¹ the calculation of ΔH_vap,G4 at 298.15 K follows eqn (3):

ΔH_vap,G4 = E_gas(T) − E_liq(T) + p(V_gas − V_liq) (2)

Given that V_liq ≪ V_gas and the kinetic energies t_gas = t_liq according to the equipartition theorem, this expression can be reduced to:

ΔH_vap,G4 = U_gas(T) − U_liq(T) + RT + C (3)

where the correction term related to vibrational motion C is very small, hence omitted in this work. The distances between CH₃ end groups were also calculated over 12 ns of this same NVT trajectory.

In parallel, an NVT simulation of 2 ns was launched using the NPT equilibrated simulation box, with a timestep of 0.5 fs. The IR spectra were generated automatically for the last 500 ps using the register_spectra keyword, from the autocorrelation of the box dipole moment. The spectra were calculated on the fly from the average of 2.0 ps segments (tseg 2.0) which corresponds to a precision of 5 cm⁻¹ on the final spectra. The influence of friction coefficient coupling the system and the thermostat was offset via the IR_Deconvolution keyword.⁸² Note that the spectra were not modified when a total time of 500 ps, 2 ns or 5 ns was chosen for the calculation. The C–H bond lengths were evaluated over the last 1 ns.

For purely classical calculations (valence fit and valence + van der Waals fit) the Berendsen barostat and the Bussi thermostat were chosen, as well as the velocity-Verlet integrator. For IR spectra exceptionally, the Langevin thermostat was used. For the semi-classical simulations, the adQTB thermostat was used in conjunction with a Langevin barostat. The BAOABRESPA integrator was employed for these simulations.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available from the authors upon reasonable request.

Supplementary information: TFSI⁻ force field parameters. Simulation results and graph analysis. See DOI: https://doi.org/10.1039/d5cp01631b

Acknowledgements

The authors thank Marin Vatin for adapting his code to systems with multiple interacting atom pairs, and Mathieu Salanne and Jean-François Dufrêche for fruitful discussions. This project was supported by the French National Research Agency under the STORE-EX Labex project (ANR-10LABX-0076) and the France 2030 program (Grant ANR-22-PEBA-0002) from the PEPR BATMAN. Calculations were performed on the HPC resources of IDRIS under the project 2023-103371.

References

1. S.-J. Tan, X.-X. Zeng, Q. Ma, X.-W. Wu and Y.-G. Guo, Electrochem. Energy Rev., 2018, 1, 113–138.
2. Z. Zhang, X. Wang, X. Li, J. Zhao, G. Liu, W. Yu, X. Dong and J. Wang, Mater. Today Sustainability, 2023, 21, 100316.
3. S. Yuan, K. Ding, X. Zeng, D. Bin, Y. Zhang, P. Dong and Y. Wang, Adv. Mater., 2023, 35, 2206228.
4. L. Elbinger, M. Enke, N. Ziegenbalg, J. C. Brendel and U. S. Schubert, Energy Storage Mater., 2024, 65, 103063.
5. A. K. X. Tan and S. Paul, Energies, 2024, 17, 5768.
6. J. W. Choi and D. Aurbach, Nat. Rev. Mater., 2016, 1, 16013.
7. T. Placke, R. Kloepsch, S. Dühnen and M. Winter, J. Solid State Electrochem., 2017, 21, 1939–1964.
8. A. Ponrouch, J. Bitenc, R. Dominko, N. Lindahl, P. Johansson and M. Palacin, Energy Storage Mater., 2019, 20, 253–262.
9. Q. Wei, L. Zhang, X. Sun and T. L. Liu, Chem. Sci., 2022, 13, 5797–5812.
10. X. Wang, Y. Huang, L. Li, L. Huang, X. Chen and Z. Yang, Chem. Phys., 2021, 542, 111051.
11. L. H. B. Nguyen, T. Picard, C. Iojoiu, F. Alloin, N. Sergent, M.-L. Doublet and J.-S. Filhol, Phys. Chem. Chem. Phys., 2022, 24, 21601–21611.
12. A. Kopac Lautar, J. Bitenc, T. Rejec, R. Dominko, J.-S. Filhol and M.-L. Doublet, J. Am. Chem. Soc., 2020, 142, 5146–5153.
13. G. R. Pastel, T. P. Pollard, O. Borodin and M. A. Schroeder, Chem. Rev., 2025, 125, 3059–3164.
14. P. Kubisiak and A. Eilmes, J. Phys. Chem. B, 2020, 124, 9680–9689.
15. P. Wróbel, P. Kubisiak and A. Eilmes, J. Phys. Chem. B, 2021, 125, 10293–10303.
16. K.-I. Saitoh, Y. Takai, T. Sato, M. Takuma and Y. Takahashi, Batteries, 2022, 8, 27.
17. Y. Zhang, E. Carino, N. T. Hahn, N. Becknell, J. Mars, K. S. Han, K. T. Mueller, M. Toney, E. J. Maginn and S. Tepavcevic, J. Phys. Chem. Lett., 2023, 14, 11393–11399.
18. N. J. Shah, C. Fang, N. C. Osti, E. Mamontov, X. Yu, J. Lee, H. Watanabe, R. Wang and N. P. Balsara, Nat. Mater., 2024, 23, 664–669.
19. Z. Tian, Y. Zou, G. Liu, Y. Wang, J. Yin, J. Ming and H. N. Alshareef, Adv. Sci., 2022, 9, 2201207.
20. M. A. Salvador, R. Maji, F. Rossella, E. Degoli, A. Ruini and R. Magri, Batteries, 2023, 9, 234.
21. L. Scalfi, M. Salanne and B. Rotenberg, Annu. Rev. Phys. Chem., 2021, 72, 189–212.
22. X. Shan, Y. Zhong, L. Zhang, Y. Zhang, X. Xia, X. Wang and J. Tu, J. Phys. Chem. C, 2021, 125, 19060–19080.
23. F. Chen, Curr. Opin. Electrochem., 2022, 35, 101086.
24. C. Zhang, J. Cheng, Y. Chen, M. K. Y. Chan, Q. Cai, R. P. Carvalho, C. F. N. Marchiori, D. Brandell, C. M. Araujo, M. Chen, X. Ji, G. Feng, K. Goloviznina, A. Serva, M. Salanne, T. Mandai, T. Hosaka, M. Alhanash, P. Johansson, Y.-Z. Qiu, H. Xiao, M. Eikerling, R. Jinnouchi, M. M. Melander, G. Kastlunger, A. Bouzid, A. Pasquarello, S.-J. Shin, M. M. Kim, H. Kim, K. Schwarz and R. Sundararaman, J. Phys. Energy, 2023, 5, 041501.
25. L. Yu, X. Chen, N. Yao, Y.-C. Gao and Q. Zhang, J. Mater. Chem. A, 2023, 11, 11078–11088.
26. E. Jónsson, Energy Storage Mater., 2020, 25, 827–835.
27. N. Yao, X. Chen, Z.-H. Fu and Q. Zhang, Chem. Rev., 2022, 122, 10970–11021.
28. N. N. Rajput, T. J. Seguin, B. M. Wood, X. Qu and K. A. Persson, in Modeling Electrochemical Energy Storage at the Atomic Scale, Springer International Publishing, Cham, 2018, pp. 79–124.
29. G. Åvall, J. Mindemark, D. Brandell and P. Johansson, Adv. Energy Mater., 2018, 8, 1703036.
30. B. Boz, T. Dev, A. Salvadori and J. L. Schaefer, J. Electrochem. Soc., 2021, 168, 090501.
31. J. N. Canongia Lopes and A. A. Pádua, Theor. Chem. Acc., 2012, 131, 1–11.
32. D. Bedrov, J.-P. Piquemal, O. Borodin, A. D. MacKerell, B. Roux and C. Schröder, Chem. Rev., 2019, 119, 7940–7995.
33. R. Fuentes-Azcatl and M. González-Melchor, J. Mol. Liq., 2023, 391, 123343.
34. X. He, Y. Zhang, H. Li, M. Zhou, W. Wang, R. Wang, K. Jiang and K. Wang, J. Energy Chem., 2025, 101, 34–41.
35. M. V. Fedkin, Y. K. Shin, N. Dasgupta, J. Yeon, W. Zhang, D. van Duin, A. C. T. van Duin, K. Mori, A. Fujiwara, M. Machida, H. Nakamura and M. Okumura, J. Phys. Chem. A, 2019, 123, 2125–2141.
36. S. Gong, Y. Zhang, Z. Mu, Z. Pu, H. Wang, Z. Yu, M. Chen, T. Zheng, Z. Wang and L. Chen, et al., arXiv, 2024, preprint, arXiv:2404.07181 DOI:10.48550/arXiv:2404.07181.
37. X. He, B. Walker, V. H. Man, P. Ren and J. Wang, Curr. Opin. Struct. Biol., 2022, 72, 187–193.
38. A. Z. Panagiotopoulos and S. Yue, J. Phys. Chem. B, 2023, 127, 430–437.
39. Z. Luo, S. A. Burrows, S. K. Smoukov, X. Fan and E. S. Boek, J. Phys. Chem. B, 2023, 127, 2224–2236.
40. J. G. McDaniel, E. Choi, C. Y. Son, J. R. Schmidt and A. Yethiraj, J. Phys. Chem. B, 2016, 120, 7024–7036.
41. J. G. McDaniel and A. Yethiraj, J. Phys. Chem. B, 2019, 123, 3499–3512.
42. K. Goloviznina, J. N. Canongia Lopes, M. Costa Gomes and A. A. H. Pádua, J. Chem. Theory Comput., 2019, 15, 5858–5871.
43. K. Goloviznina, Z. Gong, M. F. Costa Gomes and A. A. H. Pádua, J. Chem. Theory Comput., 2021, 17, 1606–1617.
44. N. Wang and E. J. Maginn, J. Phys. Chem. B, 2024, 128, 871–881.
45. M. Rezaei, S. Sakong and A. Groß, J. Chem. Theory Comput., 2023, 19, 5712–5730.
46. O. Borodin, J. Phys. Chem. B, 2009, 113, 11463–11478.
47. J. M. Litman, C. Liu and P. Ren, J. Chem. Inf. Model., 2022, 62, 79–87.
48. E. A. Vázquez-Montelongo, J. E. Vázquez-Cervantes and G. A. Cisneros, Int. J. Mol. Sci., 2020, 21, 697.
49. A. Pierini, V. Piacentini, J. L. Gómez-Urbano, A. Balducci, S. Brutti and E. Bodo, ChemistryOpen, 2024, 13, e202400134.
50. O. N. Starovoytov, J. Phys. Chem. B, 2021, 125, 11242–11255.
51. T. Binninger, D. Saraç, L. Marsh, T. Picard, M.-L. Doublet and C. Raynaud, J. Chem. Theory Comput., 2023, 19, 1023–1034.
52. T. Binninger, D. Saraç, L. Marsh, T. Picard, M.-L. Doublet and C. Raynaud, J. Chem. Theory Comput., 2025, 21, 1015.
53. T. E. Markland and M. Ceriotti, Nat. Rev. Chem., 2018, 2, 0109.
54. M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie, A. Michaelides, M. A. Morales and T. E. Markland, Chem. Rev., 2016, 116, 7529–7550.
55. N. Stolte, J. Daru, H. Forbert, J. Behler and D. Marx, J. Phys. Chem. Lett., 2024, 15, 12144–12150.
56. O. Marsalek and T. E. Markland, J. Phys. Chem. Lett., 2017, 8, 1545–1551.
57. G. R. Medders and F. Paesani, J. Chem. Theory Comput., 2015, 11, 1145–1154.
58. S. Habershon, D. E. Manolopoulos, T. E. Markland and T. F. Miller, Annu. Rev. Phys. Chem., 2013, 64, 387–413.
59. H. Dammak, Y. Chalopin, M. Laroche, M. Hayoun and J.-J. Greffet, Phys. Rev. Lett., 2009, 103, 190601.
60. K. Zhou and B. Liu, Molecular Dynamics Simulation, Elsevier, 2022, pp. 67–96.
61. F. Brieuc, Y. Bronstein, H. Dammak, P. Depondt, F. Finocchi and M. Hayoun, J. Chem. Theory Comput., 2016, 12, 5688–5697.
62. E. Mangaud, S. Huppert, T. Plé, P. Depondt, S. Bonella and F. Finocchi, J. Chem. Theory Comput., 2019, 15, 2863–2880.
63. National Institute of Advanced Industrial Science and Technology, SDBSWeb, https://sdbs.db.aist.go.jp, accessed: 29-12-2023.
64. D. Kodama, M. Kanakubo, M. Kokubo, S. Hashimoto, H. Nanjo and M. Kato, Fluid Phase Equilib., 2011, 302, 103–108.
65. C. Riadigos, R. Iglesias, M. Rivas and T. Iglesias, J. Chem. Thermodynamics, 2011, 43, 275–283.
66. G. Nichols, J. Orf, S. M. Reiter, J. Chickos and G. W. Gokel, Thermochim. Acta, 2000, 346, 15–28.
67. S. H. Jamali, T. V. Westen, O. A. Moultos and T. J. H. Vlugt, J. Chem. Theory Comput., 2018, 14, 6690–6700.
68. Y. Zhang, Y. Zhang, M. J. McCready and E. J. Maginn, J. Chem. Eng. Data, 2018, 63, 3488–3502.
69. Y. Qiu, W. Shan and H. Zhang, J. Chem. Inf. Model., 2021, 61, 3571–3582.
70. N. T. Hahn, J. Self, D. M. Driscoll, N. Dandu, K. S. Han, V. Murugesan, K. T. Mueller, L. A. Curtiss, M. Balasubramanian, K. A. Persson and K. R. Zavadil, Phys. Chem. Chem. Phys., 2022, 24, 674–686.
71. M. Vatin, M. Duvail, P. Guilbaud and J.-F. Dufrêche, J. Phys. Chem. B, 2021, 125, 3409–3418.
72. N. Mauger, T. Plé, L. Lagardère, S. Huppert and J.-P. Piquemal, J. Phys. Chem. B, 2022, 126, 8813–8826.
73. S. Bi and M. Salanne, Chem. Sci., 2024, 15, 10908–10917.
74. N. N. Rajput, X. Qu, N. Sa, A. K. Burrell and K. A. Persson, J. Am. Chem. Soc., 2015, 137, 3411–3420.
75. N. Molinari, J. P. Mailoa and B. Kozinsky, Chem. Mater., 2018, 30, 6298–6306.
76. V. Prabhakaran, G. Agarwal, J. D. Howard, S. Wi, V. Shutthanandan, D.-T. Nguyen, L. Soule, G. E. Johnson, Y.-S. Liu, F. Yang, X. Feng, J. Guo, K. Hankins, L. A. Curtiss, K. T. Mueller, R. S. Assary and V. Murugesan, ACS Appl. Mater. Interfaces, 2023, 15, 7518–7528.
77. A. C. Patel, S. Sinha and G. Palermo, Q. Rev. Biophys., 2024, 57, e15.
78. S. Bougueroua, Y. Aboulfath, A. Cimas, A. Hashemi, E. A. Pidko, D. Barth and M.-P. Gaigeot, C. R. Chim, 2024, 27, 1–23.
79. L. Lagardère, L.-H. Jolly, F. Lipparini, F. Aviat, B. Stamm, Z. F. Jing, M. Harger, H. Torabifard, G. A. Cisneros, M. J. Schnieders, N. Gresh, Y. Maday, P. Y. Ren, J. W. Ponder and J.-P. Piquemal, Chem. Sci., 2018, 9, 956–972.
80. J. A. Rackers, Z. Wang, C. Lu, M. L. Laury, L. Lagardère, M. J. Schnieders, J.-P. Piquemal, P. Ren and J. W. Ponder, J. Chem. Theory Comput., 2018, 14, 5273–5289.
81. J. Wang and T. Hou, J. Chem. Theory Comput., 2011, 7, 2151–2165.
82. M. Rossi, V. Kapil and M. Ceriotti, J. Chem. Phys., 2018, 148, 102301.

---

Footnote

† Present address: CEA, DES, ISEC, DMRC, Univ Montpellier Bagnols-sur-cèze, France.

---