Evidence for multi-Li + clustering in lithium-ion battery electrolytes

Over the past 10 years there has been intense interest in carrying out molecular dynamics simulations of lithium-ion battery (LIB) electrolytes, which in principle can yield a complete picture of electrolyte structure and dynamics. However, traditional molecular mechanics force fields, such as those used successfully for biomolecular simulation of protein–ligand binding^25,26,27, may not be well suited to describing the complexities of LIB chemistry. Specifically, the physical chemistry of Li⁺ solvation and dynamics in concentrated, polar organic electrolytes is vastly different from that of aqueous biomolecular species. Ab initio molecular dynamics (AIMD) approaches are an alternative, but suffer from severe limitations in the system sizes, simulation timescales and accuracy of the density functional theory (DFT) functional that can be accessed^28,29.

The past decade has seen the creation and evolution of a methodology for constructing force fields for complex chemical systems, which holds the promise of replicating the accuracy of high-level ab initio quantum chemistry, while reducing the cost of energy evaluation per configuration of atoms by many orders of magnitude^4,5,6,7. These methods are based on the use of machine-learning force fields (MLFFs), which are parameterized to reproduce quantum-chemical relative energies for a large and diverse training set, for example, clusters taken from force-field simulations^{7,8,9,10,11,12} or directly from periodic electronic structure calculations^13,14,15,16. The computational cost of the MLFF is two to three orders of magnitude more expensive per energy evaluation than a traditional pairwise additive molecular mechanics force field, such as optimized potentials for liquid simulations-4 (OPLS4)²⁶; however, it is still many orders of magnitude less expensive than a purely quantum-chemical approach as embodied in AIMD.

In an initial study of LIB electrolytes with a high-dimensional, charge-recursive neural network (QRNN) MLFF, molecular complexes containing multiple Li⁺ ions at short distances (as small as 3 Å) were observed without direct coordination with counterions⁹. In the present work, we report a detailed theoretical analysis of refined MLFF electrolyte model based on QRNN, coupled with comparison with solid-state nuclear magnetic resonance (ssNMR) spectroscopy. These results strongly support the computational discovery of multi-Li⁺ species exclusively coordinated by ethylene carbonate (EC) solvent molecules. We rigorously demonstrate the agreement of QRNN with high-level quantum-chemical calculations and with experimental condensed-phase, structural and electrochemical properties. Insights into the microsolvation structure of Li⁺ ions and their associated redox properties in the battery electrolyte system emerge from the QRNN molecular dynamics simulations (Fig. 1a).

Significant Li⁺ ion clustering was observed in QRNN molecular dynamics simulations of lithium hexafluorophosphate (LiPF₆) in EC at various concentrations. To quantify the liquid electrolyte structure of 1 M LiPF₆ in EC, we computed the radial distribution functions (RDFs), g(r), where r is the pair distance, of Li⁺ with other Li⁺ ions (Li–Li), with the centre of mass of EC (Li–EC), and with the centre of mass of PF₆⁻ (Li–P). These plots are shown in Fig. 1c. The Li–Li RDF from the 1 M QRNN trajectory shows a sharp peak around 3 Å, representative of the presence of multi-Li⁺ clusters in electrolyte solution, which is not seen in the OPLS4 simulations (matching results in the literature for other versions of OPLS³⁰). The Li–EC RDFs from the QRNN and OPLS trajectories are similar, representing first-shell binding of solvent to Li, but the Li–P RDFs diverge considerably, reflecting the separation of the Li⁺–PF₆⁻ ion pairs in the QRNN simulations as noted above.

The statistics pertaining to the molecular complexation of Li⁺ in different concentrations of LiPF₆ in EC, specifically the fractions of Li⁺ in different Li_n clusters are shown in Fig. 1d. In the QRNN simulations, multi-Li⁺ complexes appear frequently with more than one Li cation in close proximity (<4.5 Å), increasing in frequency with concentration. This is in contrast to the OPLS4 simulations, in which the clusters overwhelmingly contain one Li⁺ ion (more than approximately 99%); a significant proportion of clusters in the QRNN simulations contain 2 or 3 Li⁺ ions (some contain 4 or 5 Li⁺ ions, although with lower frequency, except at 2 M concentration). At 1.2 M LiPF₆, we see that 40% of the Li⁺ ions are present in multi-Li⁺ clusters during the QRNN trajectory. This translates to an average of 1.3 Li⁺ ions per Li⁺ complex in solution. Furthermore, the proportion of multi-Li⁺ complexation in EC appears to scale with concentration. At 2 M LiPF₆, more than 68% of the Li⁺ ions are in multi-Li⁺ complexes whereas at 0.5 M LiPF₆ that proportion is around 23%. The complexation statistics shown in terms of total cluster charge and categorized using the contact ion pair and aggregate nomenclature can be found in Supplementary Figs. 13 and 14. An aggregate is defined as a cluster containing more than two ions, with at least one Li⁺ and at least one PF₆⁻.

The above data clearly show that the solvation structure predicted in our recent QRNN molecular dynamics is radically different from that observed in other molecular dynamics simulations, which determines the observed battery properties. The following sections demonstrate the validity of the QRNN results via quantum-chemical calculations, macroscopic electrolyte property predictions, and comparison of experimental and theoretical NMR spectra.

Overview

Below we describe three complementary approaches to the validation of our QRNN model for the LIB electrolyte. First, we assess the accuracy with which the MLFF reproduces accurate quantum-chemical calculations for Li–Li distance scans and for the energetics of structures extracted from the QRNN and OPLS4 simulations. These results are contrasted with a similar assessment for OPLS4. Second, QRNN simulations are benchmarked for a wide range of condensed-phase properties in comparison with OPLS4, a 2022 version of QRNN, and experimental data. It is emphasized that no experimental input was used in constructing the QRNN potential energy surface, making this a fully independent test set. Finally, we demonstrate that the structure of the electrolyte from QRNN yields a ¹³C NMR signal in strong agreement with accompanying experiments across three LiPF₆ concentrations.

Ab initio validation of QRNN

First, we selected representative subsets of the Li-centred clusters from both QRNN and OPLS4 trajectories (necessary as the structures are radically different) for further validation against DFT. We calculated the gas-phase electronic energies of association for these subsets of complexes using QRNN, OPLS4 and DFT (ωB97X-D3BJ with def2-TZVPD basis set), the methodology on which the QRNN model was trained. As shown in Fig. 2b, the association energies calculated with QRNN agree within chemical accuracy (<1 kcal mol⁻¹) with those from DFT for both QRNN- and OPLS4-predicted complexes, despite their very different structures. In contrast, the errors obtained using OPLS4 are 10–100-times larger, and OPLS4 does particularly poorly at calculating the energies of the multi-Li⁺ clusters. Validation of our DFT methodology to within a few kcal mol⁻¹ against coupled cluster with singles, doubles and perturbative triples (CCSD(T)) and auxiliary-field quantum Monte Carlo (AFQMC) calculations is presented in Supplementary Information, section 2 and in ref. ³¹, respectively.

We also considered Li⁺–Li⁺ distance scans in the presence of solvent molecules. We created a test set of five such scans, three with pure EC and two with an EC/dimethyl carbonate (DMC) mixture, by randomly packing solvent molecules and ions. DMC was introduced to test whether minima would persist in the presence of a less polarizable linear carbonate. Each scan is relaxed first with OPLS then with QRNN, and finally evaluated with DFT at the ωB97X-D3BJ/def2-TZVPD level. The results are shown in Fig. 2d. A clear minimum-energy configuration emerges at a distance between 2.4 Å and 3.0 Å for each scan, which QRNN correctly reproduces. These minima, and the fact that QRNN correctly reproduces them, provide validation of the existence of low-energy Li⁺–Li⁺ complexation by solvent molecules.

Finally, we carry out an AIMD simulation on a 1.2 M LiPF₆ in EC system of nearly 500 atoms, to interrogate the stability of a di-lithium cluster in the context of fully ab initio solvation. The results show that the Li–Li distance remains in the minimum region over the duration of the simulation; details are given in Supplementary Information, section 7.

Simulation of condensed-phase properties

As in our earlier work⁹, we can validate our model by simulating a number of standard condensed-phase properties and comparing with experiment. Figure 2c shows results for predicted densities of seven pure solvents and viscosity calculations as a function of LiPF₆ concentration for a 3:7 EC/DMC mixture. Tabulated values can be found in Supplementary Tables 7 and 8. A marked improvement in accuracy and robustness compared with both the results of our previous work⁹ and those obtained from OPLS4 simulations is apparent. Our results are improved by the addition of more quantum-chemical data and better training methods in producing the QRNN model, and also by the use of a faster implementation of QRNN in the Desmond molecular dynamics engine³², which allows the use of larger simulation boxes and longer timescales. Details of the methodological improvements to both the simulation engine and QRNN development effort are given in Methods. In addition, we calculate densities, viscosities (Supplementary Information, section 10), diffusion coefficients (Supplementary Information, section 9.3) and overall fractions of Li-coordinated EC (Supplementary Information, section 11) for 0.5 M, 1.2 M and 2.0 M, compared with experimental values. Regarding transport properties, the ionic conductivities from QRNN simulations are in strong agreement at working temperatures with experimental values for 1 M LiPF₆ in EC (Supplementary Fig. 3).

NMR characterization of electrolyte

Neat EC and electrolyte samples were flash-frozen to isolate and investigate microscopic structure with ssNMR as fast exchange in the liquid electrolyte prevents resolution of different solvation environments at temperature T = 313 K (Supplementary Information, sections A.9.2 and A.9.3). First, we used direct-polarization ⁷Li and ¹H/¹⁹F → ⁷Li cross-polarization magic-angle spinning (CPMAS) ssNMR to assess Li–H and Li–F proximity. As shown in Fig. 3, direct-polarization ⁷Li ssNMR yields the highest ⁷Li signal intensity across all electrolyte compositions, consistent with detection of the full Li population. For all three samples, we see that the signal decreases when performing ¹H → ⁷Li CPMAS, which is consistent with C=O interacting with Li⁺ and protons facing away, which leads to poor polarization (see Supplementary Information for cross-polarization efficiency calculations). Further signal attenuation is observed when we perform ¹⁹F → ⁷Li CPMAS, where ⁷Li NMR is visible only after increasing the signal intensity fivefold. Both contact ion pairs and aggregates have short Li–F distances and thus should exhibit strong ⁷Li–¹⁹F dipolar couplings that enhance ¹⁹F → ⁷Li CPMAS ssNMR. The absence of this spectroscopic signature provides compelling evidence that frozen samples contain predominantly solvent-separated species, consistent with single- and multi-Li clusters. (Supplementary Figs. 23 and 24 show efforts to assign discrete ⁶Li environments).

To gain deeper insight into the local solvation environment around Li⁺, we performed ¹H → ¹³C CPMAS ssNMR (Fig. 3d,e). In neat EC, the carbonyl carbon resonates at 157.6 ppm and the ether moiety at 68.0 ppm. After the addition of LiPF₆, both peaks broaden and develop shoulders. Resolution between the parent peak and this new shoulder is more apparent for the carbonyl line, indicating that Li ions interact with the C=O group. These shoulders become more pronounced with increasing salt concentration (from 0.5 M to 1.2 M to 2.0 M), pointing to enhanced ionic aggregation, compositionally distinct solvation structures and reduced solvent availability per cation.

Next, Li⁺-centred clusters were extracted from the QRNN trajectory, relaxed with QRNN and their NMR chemical shifts computed at the PBE/pcSseg-3 level to assign EC–Li⁺ coordination^33,34,35. Figure 3d compares experimental and simulated spectra from the QRNN trajectory, and Fig. 3e shows contributions from different structural components that contribute to the peak shape (for example, free solvent, single-Li and multi-Li clusters). We see good agreement between experiment and QRNN for the spread in chemical shift of EC molecules in the Li clusters, as well as the gradual growth of a shoulder in the carbonyl region with increasing salt concentration.

When using an identical protocol to calculate NMR spectra from OPLS4, we see that simulation can adequately predict spectra for samples that predominantly contain single Li⁺–EC structures (0.5 M), but deviate from experiment when multi-Li⁺ clusters are present, especially at carbon sites closest to Li⁺ (Fig. 3d). Owing to higher irregularity in Li⁺ coordination sites and the presence of PF₆⁻-containing clusters, OPLS4 simulations fail to reproduce the line shape of the carbonyl peak of the Li⁺-centred clusters at 1.2 M and 2.0 M. At 2.0 M, the spectrum reflects a more heterogeneous and strongly perturbed local structure, leading to increased deviation and broadening compared with the free solvent, as opposed to a free-solvent-like environment, owing to the reduced ether-oxygen participation in direct Li⁺ binding. Further insights from simulated NMR spectra are provided in Supplementary Information, section 9.

The reduction potentials of species in the LIB electrolyte are crucial to the composition and functionality of the resulting anode solid electrolyte interphase (SEI). We analysed a series of Li–EC complexes, many of which were highly statistically relevant in the electrolyte from QRNN molecular dynamics, as described above, to determine their reduction potentials versus the Li/Li⁺ redox couple via quantum-chemical calculations. Geometry optimized structures (representing the lowest-energy configurations of the oxidized and reduced species) in implicit solvent using the solvation model based on density (SMD)³⁶, at the ωB97X-D3BJ/def2-TZVPD level of theory, are utilized to compute the relevant free energies. Voltages are referenced to Li/Li⁺ using a reference FeCp⁺/0 couple. Results for a series of species observed with significant frequency in the simulations are presented in Fig. 4. Further methodological details, benchmarks of the methods for computing referenced reduction potentials and reduction potentials of additional species are presented in Supplementary Information, section 12.

From Fig. 4, it is evident that the most common Li₁ complex, LiEC₄, is reduced at a much lower voltage (less favourably) than the group of Li₂EC₄₋₇ complexes, the most statistically relevant of which is Li₂EC₇. Experimentally, SEI formation onset is observed at approximately 0.8 V versus Li/Li⁺ for EC-based electrolytes on both graphite and silicon anodes, according to impedance spectroscopy and electrochemical cycling measurements^{3,17,18,19,20}. These data imply that single-Li clusters present in the electrolyte would not be spontaneously reduced at this voltage. In contrast, the redox potentials for Li₂ clusters fall within 0.1 V of this value. In addition, the Li₂ complexes appear to show even more favourable second one-electron reductions (Supplementary Table 13), adding the possibility that these complexes undergo spontaneous further reduction following the initial one-electron process around 0.8 V versus Li/Li⁺. Altogether, the more favourable one-electron reduction potentials for multi-Li⁺ complexes suggest that they participate in, if not govern, the initial anode SEI formation.

Having established the presence of multi-Li⁺ clusters at typical salt concentrations used in LIBs (approximately 1 M), we now consider what role these clusters play in the operation of the cell. After polarization of the cell, we expect that positively charged species (such as those contained in the multi-Li⁺ clusters) are attracted to the negatively charged anode surface in the electrochemical double layer. These species are probably the reactants that reduce at 0.8 V, forming the first layer of the SEI, and are responsible for the products reported in literature (lithium ethylene dicarbonate, Li₂CO₃, lithium butylene dicarbonate, oligomers and so on). We also speculate that the multi-Li⁺ clusters are involved in the transport of Li ions from cathode to anode and the release of these ions at the electrode/electrolyte interface during battery charging.

SEI formation is crucial to proper LIB functioning, and from the electrochemical calculations above, it is reasonable to hypothesize that the reductive decomposition reactions necessary to SEI formation proceed through the multi-Li⁺ cluster species. This deduction in turn provides an explanation as to why EC and related solvents are utilized in the LIB electrolyte. It is almost certainly the case that the high dielectric constant and polarizability of EC (unusual in an aprotic solvent) have major roles in stabilizing the co-complexation of multiple Li⁺ ions. In essence the EC molecules coordinated to the Li⁺ ions of the cluster provide a potent microsolvation, enabling close approach of the positive ions, which would be impossible in the gas phase or in a normal organic solvent. The bridging location of the EC carbonyl oxygen is probably crucial. Further analysis of the decomposition of the multi-Li clusters in the presence of the electrode via QRNN simulation offers the promise of achieving a detailed atomic-level picture of the formation of the SEI and its dependence on electrolyte composition. This capability in turn opens the door to using large-scale graphics-processing-unit-based simulation to carry out electrolyte engineering, thus enabling the design of more efficient and cost-effective LIBs and, in principle, batteries made from other constituents as well.

The above analyses of the LiPF₆ in EC electrolyte represent a substantial step forward in understanding LIB functionality at an atomic level of detail, providing surprising and previously unknown insights into the microscopic structure of the electrolyte. In future work, we intend to use the QRNN simulation platform to study a wide range of electrolyte solvents, salts and additives, comparing with experimental data when available. We will furthermore add electrode models, a description of reactive chemistry and a treatment of the electrode applied voltage, enabling a high degree of realism in the investigation of the early stages of SEI formation to be attained.

In addition to specific implications for understanding the structure and function of EC-based LIB electrolytes, the present work provides proof of concept that quite complex and even surprising phenomena can be accurately modelled via a properly parameterized MLFF without any explicit experimental input. Previous work on water has shown similar capabilities^37,38, but the LIB electrolyte system is far more complex in composition and in the forming of unusual covalently bonded species as opposed to the standard non-bonded interactions (van der Waals and electrostatic) that dominate liquid water. The results of this paper suggest that a wide range of important chemical systems, composed of elements from across the full range of the periodic table, can now be effectively addressed via computationally accessible dynamical molecular simulations.

MLFF and molecular dynamics

Our simulations were performed at 313 K in the canonical (NVT) ensemble (in which the number of particles (N), volume (V) and temperature (T) are kept constant), above room temperature because pure EC is a solid at room temperature. For QRNN, the systems were equilibrated for 1 ns in an isothermal–isobaric (NPT) ensemble (in which the number of particles (N), pressure (P) and temperature (T) are kept constant) before 5 ns of production, whereas for the OPLS4 molecular dynamics, we performed an equilibration in NPT for 10 ns followed by 50 ns of production dynamics in NVT. The Li-based molecular clusters were extracted from molecular dynamics frames using a 4.5-Å radius around Li, including molecules with any atoms within this radius.

We carried out molecular dynamics simulations of 1 M LiPF₆ solution in EC solvent using both a machine-learned force field (QRNN) and an empirical force field (OPLS4). We ran NVT trajectories at 313 K for 5 ns using the refined QRNN force field and 50 ns using the OPLS4 force field. Molecular complexes were extracted from these trajectories at evenly spaced time steps based on the immediate coordination environment of the Li ions in each simulation. The equilibrium dynamics of the same electrolyte yielded vastly different molecular structures depending on the force field used leading to previously unforeseen compositions of Li complexes.

A batch of replicas is launched from the command line onto Google Cloud resources by specifying the desired system composition. Each replica is automatically assigned a unique random seed, which is used to construct an initial simulation cell of approximately 6,000 atoms using a Monte Carlo packing algorithm implemented in Schrödinger’s Disordered System Builder³². This ensures that each replica has an independent starting point for maximum sampling efficiency. The replicas proceed automatically to an equilibration stage consisting of 1 ns (2 million steps) of NPT molecular dynamics with our Desmond QRNN molecular dynamics engine³².

After equilibration, instead of simply using the final simulation frame for the following steps, our scripts automatically select the trajectory frame with a volume closest to the mean volume over the final 30% of each trajectory. This ensures that the density used for the following step is as close to the true mean density of the model system as possible. After equilibration, the replicas proceed to 5 ns (10 million steps) of NVT molecular dynamics for property calculation.

As in our previous MLFF papers^5,6,9,10,11, our model was trained using a fork of the torchANI library^40,41, with data labelled at the ωB97X-D3BJ/def2-TZVPD level of theory^31,42,43 using Psi4 (ref. ⁴⁴) with the same settings as in ref. ¹⁰:

'scf_type': 'MEM_DF', 'dft_basis_tolerance': 1 × 10⁻¹⁰

'ints_tolerance': 1e-10, 'maxiter': 200, 'dft_pruning_scheme': 'robust'

's_orthogonalization': 'partialcholesky', 's_cholesky_tolerance': 1 × 10⁻⁶

Our previous data from ref. ⁹ were checked using these settings, and inconsistent self-consistent-field solutions were removed. As in ref. ⁹, we trained the charge model using DFT dipole moments, rather than xTB atomic partial charges⁴⁵. Small-scale geometry optimization and molecular dynamics for data sampling were performed using the Atomic Simulation Environment library⁴⁶ via the Schrödinger Suite³².

Sample preparation and NMR experiments

All samples, including 0.5 M, 1.2 M and 2.0 M LiPF₆ in EC, as well as pure EC solvent, were prepared and packed into 2.5-mm Phoenix zirconia rotors entirely inside an argon-filled glovebox to maintain an inert environment and prevent exposure to air and moisture. Samples were flash-frozen before insertion into the probe to preserve their room-temperature structure during low-temperature measurements.

¹³C ssNMR experiments were performed on a Bruker Ascend 600-MHz wide-bore compact magnet equipped with a 2.5-mm HFXY LT-DNP Phoenix probe, tuned to ¹H¹³C mode. The MAS rate was maintained at 20 kHz, and the sample temperature during measurements was estimated to be approximately 116 K. The 90° pulse powers were 172 kHz at 60 W for ¹H and 75 kHz at 60 W for ¹³C. One-dimensional ¹³C spectra were recorded using a standard ¹H–¹³C cross-polarization experiment with a cross-polarization contact time of 3 ms and recycle delay of 300 s. Adamantane was used as an external chemical-shift reference for ¹³C. The chemical shift of the CH₂ peak (37.75 ppm) was adjusted for 116 K using the temperature correction formula reported in ref. ⁴⁷. Bulk ¹H T₁ relaxation times were measured using a saturation recovery experiment with cross-polarization detection on ¹³C. T₁ values were determined to be 6.3 s for the 2.0 M LiPF₆ samples and 42.2 s for the 1.2 M LiPF₆ samples. For the 0.5 M LiPF₆ and pure EC samples, ¹H T₁ could not be accurately determined as the signal did not saturate even after a 500-s delay.

The differences in ¹H T₁ between the samples may impact the relative intensities of the free solvent peaks and the solvent–ion coordination peaks in ¹H–¹³C CPMAS ssNMR. For example, we see a broad shoulder appear in ¹³C ssNMR after the addition of 0.5 M LiPF₆ to EC, indicating that Li ions coordinate to the solvent. In principle, the integrated area of this peak indicates the total amount of solvent coordinated to Li. However, the intensity of this shoulder changes very little after the salt concentration is increased to 1.2 M LiPF₆, where we expect to have additional salt–solvent interactions, probably due to differences in ¹H T₁ and cross-polarization transfer efficiency. Direct polarization ⁷Li and ¹H → ⁷Li and ¹⁹F → ⁷Li cross-polarization MAS ssNMR measurements were performed on a Bruker Avance NEO 400-MHz wide-bore spectrometer equipped with a 2.5-mm HFXY LT-DNP Phoenix probe tuned to ¹H/¹⁹F/⁷Li mode. LiPF₆ electrolyte samples in EC at 0.5 M, 1.2 M and 2.0 M were prepared alongside an additional sample of 1.2 M LiPF₆ in EC/propylene carbonate (80:20 w/w). The MAS rate was maintained at 20 kHz and the sample temperature during measurements was approximately 210 K. The 90° pulse radio-frequency field strengths were 100 kHz (116 W) for ¹H, 75 kHz (95 W) for ¹⁹F, and 125 kHz (80 W) for ⁷Li. One-dimensional ⁷Li spectra were acquired using a single-pulse direct polarization experiment with a recycle delay of 300 s. ¹H–⁷Li cross-polarization experiments were performed with a contact time of 3 ms and a recycle delay of 10 s, and ¹⁹F–⁷Li cross-polarization experiments were performed with a contact time of 2 ms and a recycle delay of 3 s. LiF was used as an external chemical-shift reference for both ⁷Li and ¹⁹F and set to −1 ppm and −204 ppm, respectively. ⁷Li longitudinal relaxation times T₁ were measured using a saturation recovery pulse sequence and were found to be 10 s for 0.5 M, 19 s for 1.2 M, and 43 s for 2 M electrolytes (Supplementary Fig. 25). For the ⁶Li ssNMR experiments, the same samples described above were used. Experiments were performed on a Bruker Avance NEO 900-MHz spectrometer equipped with a 2.5-mm Black Fox HXY probe tuned to ¹H⁶Li mode. The MAS rate was maintained at 20 kHz and the sample temperature was set to 210 K, which was limited by the probehead available at high field that also supports ⁶Li tuning. The 90° pulse radio-frequency field strengths were 64 kHz (150 W) for ¹H and 50 kHz (150 W) for ⁶Li. One-dimensional ⁶Li spectra were acquired using a single-pulse direct polarization experiment with a recycle delay of 300 s. ¹H–⁶Li cross-polarization experiments were performed with a contact time of 8 ms and a recycle delay of 5 s. Two-dimensional ⁶Li–⁶Li radio-frequency-driven recoupling experiments were performed at 10-kHz MAS with a mixing time of 5 ms and a recycle delay of 5 s. The temperature readouts on the 2.5-mm PhoenixNMR LTMAS DNP probeheads operating at both 600 MHz and 400 MHz were calibrated using the T₁ of KBr following the procedure reported in ref. ⁴⁸. The temperature readout on the 2.5-mm Black Fox NMR probehead operating at 900 MHz was calibrated by monitoring the ¹H chemical shift from the H₆ resonance of the high-field NMR thermometer, thulium 1,4,7,10-tetra-azacyclododecane-1,4,7,10-tetrakis (methylene phosphonate) (TmDOTP_-5, often denoted TmDOTP), as a function of temperature under MAS conditions, following the procedure reported in ref. ⁴⁹.

Solution NMR samples were prepared in an argon-filled glovebox by heating samples to 40 °C before filling a 3.65-mm-outer-diameter fluorinated ethylene propylene tube with 300 μl of sample. This tube was snugly placed into a 5-mm-outer-diameter glass J-young tube containing roughly 100 μl of DMSO-d₆ used for locking.

¹H, ⁷Li, ¹⁷O and ¹⁹F NMR spectra were collected on a Bruker Avance III 400-MHz spectrometer using a triple resonance broadband observe probehead. Samples were kept at 40 °C inside the NMR probe using a calibrated variable temperature unit. For 90° pulse spectra, ¹H experiments used 32 scans with a recycle delay of 1 s and were referenced to tetramethylsilane at 0 ppm. ⁷Li experiments used 32 scans with a recycle delay of 5 × T₁ as determined through an inversion recovery experiment and were referenced to LiCl at 0 ppm. ¹⁷O experiments used 8,192 scans with a recycle delay of 0.01 s and were referenced to water at 0 ppm, and ¹⁹F experiments used 32 scans with a recycle delay of 1 s and were internally referenced to the PF₆⁻ doublet at −73 ppm.

For all pulse-field-gradient experiments, 16 linearly spaced gradient values between 2% and 98% were measured over 4 scans with recycle delays of 10 s, 15 s and 5 × T₁ for ¹H, ¹⁹F and ⁷Li respectively.

Viscosity measurements

Viscosity was determined using a Cannon–Fenske Routine size 75 viscometer (0.008 centistokes per second constant) to measure kinematic viscosity. The viscometer was thoroughly washed and dried in a 60 °C oven overnight and then transferred to an argon-filled glovebox, where the measurements were taken. The viscometer was charged with liquid electrolyte and then placed in a 313-K oil bath for at least 15 min before kinematic viscosity measurements. Efflux time measurements were repeated at least three times for each liquid sample.

To measure density, liquid samples were brought to 313 K using an oil bath. Five millilitres of each sample was massed at least five times to determine the density. Viscosity was calculated as the product of measured kinematic viscosity and density and used to calculate diffusivity from the pulse-field-gradient experiments.

Raman spectroscopy measurements

Raman spectra were collected for LiPF₆ dissolved in EC at concentrations of 0.5 M, 1.2 M and 2.0 M. Electrolyte solutions were heated to 313 K to ensure complete dissolution before transfer into glass capillary tubes in an argon-filled glovebox. The capillaries were flame-sealed at one end and sealed at the other using a malleable sealing compound to prevent air and moisture exposure. The filled capillaries were then mounted on the stage of a Renishaw inVia Raman spectrometer and aligned under the microscope objective. Spectra were acquired at ambient temperatures using a 633-nm excitation laser with an 1,800 l mm⁻¹ grating over a spectral window of 440–1,000 cm⁻¹, with an acquisition time of 3 s and 30 accumulations per spectrum at 100% laser power. The instrument was calibrated before measurement using a silicon reference standard (520 cm⁻¹ peak).

NMR simulation details

To assign EC–Li⁺ coordination structures detected in ¹³C ssNMR, NMR calculations were carried out by extracting all Li⁺-centred clusters (here called coordinated Li or Li⁺-centred clusters) from the QRNN trajectory, followed by relaxing the cluster geometries with QRNN and computing NMR chemical shifts using PBE/pcSseg-3 (refs. ^33,34,35). For all salt concentrations (0.5 M, 1.2 M and 2.0 M), the predicted NMR spectrum, excluding the free-solvent contribution, was obtained by summing individual simulated spectra for all Li⁺-centred clusters in the molecular dynamics simulation box⁵⁰. The free-solvent contribution was modelled as a single Gaussian line shape. The linewidths of the individual Li⁺-centred cluster peak, and that of the free-solvent peak, are fit to reproduce the experimental data.

To compare the full simulated NMR line shapes with experiment, it is necessary to place the free-solvent peak relative to the Li⁺-centred cluster peak. For the 0.5 M and 1.2 M concentrations, for both the carbonyl and ether peaks, the relative free-solvent–Li⁺-centred cluster peak separation is directly obtained from the experimental spectra. For the 2.0 M system, we do not average absolute chemical shifts; instead, we approximate the free-solvent–Li⁺-centred cluster peak separation by averaging the corresponding values from the 0.5 M and 1.2 M systems, which are very similar. This approximation is motivated by the fact that the free-solvent peak is less well defined at 2.0 M and that a first-principles determination of its position is not feasible, owing to the lack of a well-characterized structural model for the free-solvent region at low temperature. The approach implicitly assumes that the free-solvent–Li⁺-centred cluster peak separation varies weakly with concentration and that the dominant local environments contributing to the free-solvent region are not markedly altered between 1.2 M and 2.0 M. On the basis of the small variation between the 0.5 M and 1.2 M values, we estimate an uncertainty of approximately the same order in the free-solvent–Li⁺-centred cluster separation at 2.0 M. Importantly, the simulated line shapes and resulting spectral assignments are not sensitive to this uncertainty. As further support, in Supplementary Information, section 9.7, we present a physically motivated model for the free-solvent region that yields calculated spectra very similar to those obtained using the empirically derived free-solvent–Li⁺-centred cluster separation. Details of the simulations are discussed further in Supplementary Information, section 9.

Comparisons between the experimental and simulated spectra from the QRNN trajectory are shown in Fig. 3d; contributions from different structural components that comprise peak shape, such as free-solvent EC, single-Li and multi-Li clusters in the QRNN NMR simulation are shown in Fig. 3e. The spread in chemical shift of EC molecules in the Li clusters, as well as the gradual growth of the shoulder of the carbonyl peak with increasing salt concentration, show good agreement with experiments with the QRNN spectra. Moreover, using the average gap obtained from the 0.5 M and 1.2 M concentrations yields great agreement of the QRNN spectra with experiment across carbonyl and ether peaks for the 2 M concentration.

All the NMR calculations have been performed using ORCA 6.0 (ref. ⁵¹). The NMR chemical shielding values are calculated using PBE/pcSseg-3 (refs. ^33,34). The isotropic chemical shifts obtained are used to simulate the spectra of different cluster compositions.

Two distinct types of species are identified in the molecular dynamics trajectories, namely, coordinated Li clusters and free-solvent molecules. We follow the following steps to simulate the NMR spectra at different salt concentrations using the QRNN or OPLS4 trajectories.

(1) Coordinated Li clusters are extracted from the molecular dynamics trajectories. For the purpose of NMR calculations, the effect of PF₆⁻ ions, even when they are not directly part of the Li⁺ complexation environment but are present nearby, are considered. To achieve this, when applicable, one PF₆⁻ unit is included beyond the regular cluster when a distance cut-off is satisfied. When there are multiple qualifying PF₆⁻ ions, the one closest to the original cluster is included. NMR spectra are computed for each cluster by summing Lorentzian line-shape contributions arising from the isotropic chemical shifts of all carbon atoms within the clusters.

(2) The net contribution of all the coordinated Li clusters in the molecular dynamics trajectory is obtained by summing the spectra of each individual Li cluster.

(3) The number of EC molecules within coordinated Li clusters is counted along with the number of free-solvent EC molecules in the molecular dynamics trajectory.

(4) For the 0.5 M and 1.2 M salt concentrations, the experimental spectra is fitted by a Voigt profile each for the solvent peak and the coordinated Li peak⁵². The relative chemical shift between the two peaks is extracted.

(5) The gap between the peaks is used to calculate the position of the free-solvent peak with respect to the coordinated Li peak, calculated earlier, which is modelled by a single Gaussian line shape.

(6) The width of the free-solvent peak, as well as the individual Lorentzian line shape, is tuned to fit the experimental NMR spectra for different peaks at different concentrations.

(7) The reference is chosen so that the free solvent peak is aligned with the experimental solvent peak. For each concentration, the reference for carbonyl and ether peak is chosen individually⁵³.

(8) The sum of the coordinated Li peak and the free-solvent peak produces the final simulated spectra, to be compared with the experimental NMR spectra. The coordinated Li and free-solvent components in each simulated spectrum are intensity normalized weighted by the fraction of EC molecules in Li clusters relative to those in the free solvent in the molecular dynamics simulation box.

(9) For the 2 M concentration, instead of calculating the gap between the coordinated Li and the free-solvent peak from the experimental spectra, we use the average obtained from 0.5 M and 1.2 M and use the same gap relative to the coordinated Li peak, accounting for the relative positions of the coordinated Li peaks for 2 M concentration with respect to 1.2 M concentration.

Reduction potentials versus Li/Li⁺ via ferrocenium/ferrocene

To compute the reduction potentials of Li–electrolyte complexes on the LIB-relevant reference scale versus Li/Li⁺, we utilized experimental data and computed potentials for the ferrocenium/ferrocene redox couple (denoted hereafter as FeCp⁺/0) in LIB-relevant electrolytic conditions. We note that this approach is a departure from a considerable body of computational literature in which a constant, vacuum-based reference of 1.4 V is applied to the computed redox potentials^54,55,56,57. Our approach is, however, consistent with emerging best practices of experimental electrochemistry, in which using FeCp⁺/0 is recommended given its versatile solubility and electrochemical reversibility in many different solvents and electrolytes^57,58. Given the challenging nature of the transition-metal redox chemistry of FeCp⁺/0, we furthermore impose a correction for the error in the gas-phase adiabatic ionization potential of each electronic structure method used⁵⁹.

Thus, the reduction potentials of Li–electrolyte coordination complexes referenced to Li/Li⁺, E_red, are computed as follows:

where E_ref(FeCp^+/0) is the experimental reduction potential of FeCp^+/0 relative to Li/Li⁺ in the electrolytic conditions of interest, E_calc(FeCp^+/0) is the calculated, unreferenced reduction potential of FeCp^+/0, E is the computed reduction potential of the Li-based complex, \({\mathcal{F}}\) is Faraday's constant, n is the number of moles of electrons transferred, and δ is the error in the gas-phase adiabatic ionization potential of the electronic structure method. E_ref(FeCp^+/0) in the present work is taken to be 3.25 V versus Li/Li⁺ based on cyclic voltammetry of FeCp^+/0 in 1 M LiPF₆ in 1:1 EC/EMC⁶⁰. In this expression, E_calc(FeCp^+/0) and E are computed as the solution-phase free-energy differences between the reduced and oxidized species. As a consequence of the electrolyte conditions of the reference couple, we use ϵ = 33.6 for implicit solvation of all species in reduction potential calculations when using both conductor-like continuum polarization model (C-PCM) and SMD solvation schemes⁶¹. The additional parameters for the SMD scheme are taken as those previously determined for pure EC⁶². The values of δ for ωB97X and DLPNO-CCSD(T) are taken from previous benchmarking calculations⁵⁹ and are −4.30 kcal mol⁻¹ and +3.06 kcal mol⁻¹, respectively. The reaction free energies associated with each reduction reaction, ΔG°, were computed as the sum of the reaction free energies in gas phase, ΔG°, and the differences in SMD solvation energies between the reduced and oxidized states, ΔΔE(sol). Zero-point energy and entropic contributions to the free energies were obtained by a harmonic frequency calculation of each species at its gas-phase optimized geometry.

The total shift on the computed reduction potentials, E_ref(FeCp^+/0) + E_calc(FeCp^+/0) + δ, corresponding to the Li/Li⁺ scale are 1.36 V and 1.43 V for DLPNO-CCSD(T)/aug-cc-pVTZ and ωB97X-D3BJ/def2-TZVPD, respectively, using SMD for solvation. Although these values are close to the frequently used 1.4 V, they reflect sensitivities to the errors in the underlying methods and moreover could be consistently adopted for any electrolytic conditions as long as the accompanying experimental voltammetry for FeCp^+/0 versus Li/Li⁺ has been measured. The reduction potentials in the main text were computed using ORCA 6.0⁵¹ and the ωB97X-D3BJ/def2-TZVPD level of theory^42,43.

The spin density isosurfaces in Fig. 4 demonstrate that the excess electron associated with reduction is localized on the carbonate moiety of an EC ligand, consistent with previous DFT studies of reduced Li₁/EC complexes⁶³. In the reduced multi-Li⁺ complexes, EC ligands that coordinate two or more Li⁺ centres are the electron acceptors, which can help explain the more facile reduction of multi-Li⁺ species. The stabilization of the radical EC⁻ moiety by multiple coordinating Li⁺ cations probably contributes to a more stable reduced complex.

1. Goodenough, J. B. & Kim, Y. Challenges for rechargeable Li batteries. Chem. Mater. 22, 587–603 (2010).

   Article  CAS  Google Scholar 

2. Xu, K. Electrolytes and interphases in Li-ion batteries and beyond. Chem. Rev. 114, 11503–11618 (2014).

   Article  CAS  PubMed  Google Scholar 

3. Xu, K. Nonaqueous liquid electrolytes for lithium-based rechargeable batteries. Chem. Rev. 104, 4303–4418 (2004).

   Article  CAS  PubMed  Google Scholar 

4. Unke, O. T. et al. Machine learning force fields. Chem. Rev. 121, 10142–10186 (2021).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

5. Stevenson, J. M. et al. Schrödinger-ANI: an eight-element neural network interaction potential with greatly expanded coverage of druglike chemical space. Preprint at https://doi.org/10.48550/arxiv.1912.05079 (2019).

6. Jacobson, L. D. et al. Transferable neural network potential energy surfaces for closed-shell organic molecules: extension to ions. J. Chem. Theory Comput. 18, 2354–2366 (2022).

   Article  CAS  PubMed  Google Scholar 

7. Gong, S. et al. A predictive machine learning force-field framework for liquid electrolyte development. Nat. Mach. Intell. 7, 543–552 (2025).

8. Xie, X., Persson, K. A. & Small, D. W. Incorporating electronic information into machine learning potential energy surfaces via approaching the ground-state electronic energy as a function of atom-based electronic populations. J. Chem. Theory Comput. 16, 4256–4270 (2020).

   Article  CAS  PubMed  Google Scholar 

9. Dajnowicz, S. et al. High-dimensional neural network potential for liquid electrolyte simulations. J. Phys. Chem. B 126, 6271–6280 (2022).

10. Stevenson, J. M., Jacobson, L. D., Agarwal, G. & Dajnowicz, S. Transfer learning lithium and electrolyte potential energy surfaces from pure and hybrid DFT. In AI for Accelerated Materials Design NeurIPS 2022 Workshop (2022).

11. Stevenson, J. M., Jacobson, L. D. & Agarwal, G. Machine learning force field ranking of candidate solid electrolyte interphase structures in li-ion batteries. In AI for Accelerated Materials Design NeurIPS 2022 Workshop (2023).

12. Zeng, B. et al. Uni-ELF: A multi-level representation learning framework for electrolyte formulation design. Preprint at https://doi.org/10.48550/arxiv.2407.06152 (2024).

13. Behler, J. & Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. Phys. Rev. Lett. 98, 146401 (2007).

    Article  PubMed  ADS  Google Scholar 

14. Artrith, N., Morawietz, T. & Behler, J. High-dimensional neural-network potentials for multicomponent systems: applications to zinc oxide. Phys. Rev. B 83, 153101 (2011).

    Article  ADS  Google Scholar 

15. Zhang, L. et al. A deep potential model with long-range electrostatic interactions. J. Chem. Phys. 156, 124107 (2022).

16. Magdău, I.-B. et al. Machine learning force fields for molecular liquids: ethylene carbonate/ethyl methyl carbonate binary solvent. npj Comput. Mater. 9, 146 (2023).

    Article  Google Scholar 

17. Xu, C. et al. Improved performance of the silicon anode for li-ion batteries: understanding the surface modification mechanism of fluoroethylene carbonate as an effective electrolyte additive. Chem. Mater. 27, 2591–2599 (2015).

    Article  CAS  Google Scholar 

18. Michan, A. L., Leskes, M. & Grey, C. P. Voltage dependent solid electrolyte interphase formation in silicon electrodes: monitoring the formation of organic decomposition products. Chem. Mater. 28, 385–398 (2016).

    Article  CAS  Google Scholar 

19. Jin, Y. et al. Identifying the structural basis for the increased stability of the solid electrolyte interphase formed on silicon with the additive fluoroethylene carbonate. J. Am. Chem. Soc. 139, 14992–15004 (2017).

    Article  CAS  PubMed  ADS  Google Scholar 

20. Jin, Y. et al. Understanding fluoroethylene carbonate and vinylene carbonate based electrolytes for Si anodes in lithium ion batteries with NMR spectroscopy. J. Am. Chem. Soc. 140, 9854–9867 (2018).

    Article  CAS  PubMed  ADS  Google Scholar 

21. Väli, R., Jänes, A. & Lust, E. Vinylene carbonate as co-solvent for low-temperature mixed electrolyte based supercapacitors. J. Electrochem. Soc. 163, A851 (2016).

    Article  Google Scholar 

22. Flick, E. W. Industrial Solvents Handbook 5th edn (William Andrew Publishing, 1998).

23. Logan, E. R. et al. A study of the physical properties of Li-ion battery electrolytes containing esters. J. Electrochem. Soc. 165, A21 (2018).

    Article  CAS  Google Scholar 

24. Jänes, A., Thomberg, T., Eskusson, J. & Lust, E. Fluoroethylene carbonate and propylene carbonate mixtures based electrolytes for supercapacitors. ECS Trans. 58, 71 (2014).

    Article  ADS  Google Scholar 

25. Wang, L. et al. Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field. J. Am. Chem. Soc. 137, 2695–2703 (2015).

    Article  CAS  PubMed  ADS  Google Scholar 

26. Lu, C. et al. OPLS4: improving force field accuracy on challenging regimes of chemical space. J. Chem. Theory Comput. 17, 4291–4300 (2021).

    Article  CAS  PubMed  Google Scholar 

27. Damm, W. et al. OPLS5: addition of polarizability and improved treatment of metals. Preprint at https://doi.org/10.26434/chemrxiv-2024-2svr (2024).

28. Leung, K. & Budzien, J. L. Ab initio molecular dynamics simulations of the initial stages of solid–electrolyte interphase formation on lithium ion battery graphitic anodes. Phys. Chem. Chem. Phys. 12, 6583–6586 (2010).

    Article  CAS  PubMed  Google Scholar 

29. Soto, F. A., Marzouk, A., El-Mellouhi, F. & Balbuena, P. B. Understanding ionic diffusion through SEI components for lithium-ion and sodium-ion batteries: insights from first-principles calculations. Chem. Mater. 30, 3315–3322 (2018).

    Article  CAS  Google Scholar 

30. Hasan, R. & Datta, D. Computational study of Li⁺ solvation structures in fluorinated ether, non-fluorinated ether, and organic carbonate-based electrolytes at low and high salt concentrations. Energy Adv. 4, 1049–1066 (2025).

    Article  CAS  Google Scholar 

31. Debnath, S. et al. Accurate quantum chemical reaction energies for lithium-mediated electrolyte decomposition and evaluation of density functional approximations. J. Phys. Chem. A 127, 9178–9184 (2023).

32. Schrödinger. Schrödinger Release 2025-3: Materials Science Suite. Schrödinger https://www.schrodinger.com/platform/materials-science (2025).

33. Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 77, 3865–3868 (1996).

    Article  CAS  PubMed  ADS  Google Scholar 

34. Jensen, F. Segmented contracted basis sets optimized for nuclear magnetic shielding. J. Chem. Theory Comput. 11, 132–138 (2015).

    Article  CAS  PubMed  Google Scholar 

35. Liang, J. et al. Efficient calculation of nmr shielding constants using composite method approximations and locally dense basis sets. J. Chem. Theory Comput. 19, 514–523 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

36. Marenich, A. V., Cramer, C. J. & Truhlar, D. G. Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions. J. Phys. Chem. B 113, 6378–6396 (2009).

    Article  CAS  PubMed  Google Scholar 

37. Bore, S. L. & Paesani, F. Realistic phase diagram of water from ‘first principles’ data-driven quantum simulations. Nat. Commun. 14, 3349 (2023).

    Article  CAS  PubMed  PubMed Central  ADS  Google Scholar 

38. Chen, M. S. et al. Data-efficient machine learning potentials from transfer learning of periodic correlated electronic structure methods: liquid water at AFQMC, CCSD, and CCSD(T) accuracy. J. Chem. Theory Comput. 19, 4510–4519 (2023).

    Article  CAS  PubMed  Google Scholar 

39. Sasaki, Y. Physical and electrochemical properties and application to lithium batteries of fluorinated organic solvents. In Fluorinated Materials for Energy Conversion (eds Nakajima, T. & Groult, H.) 285–304 (Elsevier, 2005).

40. Gao, X., Ramezanghorbani, F., Isayev, O., Smith, J. S. & Roitberg, A. E. TorchANI: a free and open source PyTorch-based deep learning implementation of the ANI neural network potentials. J. Chem. Inf. Model. 60, 3408–3415 (2020).

    Article  CAS  PubMed  Google Scholar 

41. Ramezanghorbani, F. Developing Machine Learning Models to Enhance Applicability of Neural Network Potentials in Drug Discovery. PhD thesis, Univ. Florida (2020).

42. Chai, J.-D. & Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections. Phys. Chem. Chem. Phys. 10, 6615–6620 (2008).

    Article  CAS  PubMed  Google Scholar 

43. Goerigk, L. et al. A look at the density functional theory zoo with the advanced GMTKN55 database forgeneral main group thermochemistry, kinetics and noncovalent interactions. Phys. Chem. Chem. Phys. 19, 32184–32215 (2017).

    Article  CAS  PubMed  Google Scholar 

44. Turney, J. M. et al. Psi4: an open-source ab initio electronic structure program. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2, 556–565 (2012).

    Article  CAS  Google Scholar 

45. Bannwarth, C., Ehlert, S. & Grimme, S. GFN2-xTB—an accurate and broadly parametrized self-consistent tight-binding quantum chemical method with multipole electrostatics and density-dependent dispersion contributions. J. Chem. Theory Comput. 15, 1652–1671 (2019).

    Article  CAS  PubMed  Google Scholar 

46. Larsen, A. H. et al. The atomic simulation environment—a Python library for working with atoms. J. Phys. Condensed Matter 29, 273002 (2017).

    Article  Google Scholar 

47. Hoffman, R. Solid-state chemical-shift referencing with adamantane. J. Magn. Reson. 340, 107231 (2022).

    Article  CAS  PubMed  Google Scholar 

48. Thurber, K. R. & Tycko, R. Measurement of sample temperatures under magic-angle spinning from the chemical shift and spin-lattice relaxation rate of ⁷⁹Br in KBr powder. J. Magn. Reson. 196, 84–87 (2009).

    Article  CAS  PubMed  ADS  Google Scholar 

49. Zhang, D., Itin, B. & McDermott, A. E. TmDOTP: an NMR-based thermometer for magic angle spinning nmr experiments. J. Magn. Reson. 308, 106574 (2019).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

50. Atwi, R. et al. An automated framework for high-throughput predictions of NMR chemical shifts within liquid solutions. Nat. Comput. Sci. 2, 112–122 (2022).

    Article  CAS  PubMed  Google Scholar 

51. Neese, F. Software update: The orca program system—version 6.0. Wiley Interdiscip. Rev. Comput. Mol. Sci. 15, e70019 (2025).

    Article  Google Scholar 

52. Higinbotham, J. & Marshall, I. NMR lineshapes and lineshape fitting procedures. Annu. Rep. NMR Spectrosc. 43, 59–120 (2001).

53. Sarotti, A. M. & Pellegrinet, S. C. A multi-standard approach for GIAO ¹³C NMR calculations. J. Org. Chem. 74, 7254–7260 (2009).

    Article  CAS  PubMed  Google Scholar 

54. Borodin, O., Behl, W. & Jow, T. R. Oxidative stability and initial decomposition reactions of carbonate, sulfone, and alkyl phosphate-based electrolytes. J. Phys. Chem. C 117, 8661–8682 (2013).

    Article  CAS  Google Scholar 

55. Borodin, O. Molecular modeling of electrolytes. In Electrolytes for Lithium and Lithium-ion Batteries (eds. Jow, T. R. et al) 371–401 (Springer, 2014).

56. Delp, S. A. et al. Importance of reduction and oxidation stability of high voltage electrolytes and additives. Electrochim. Acta 209, 498–510 (2016).

    Article  CAS  ADS  Google Scholar 

57. Gibson, L. D., Pfaendtner, J. & Mundy, C. J. Probing the thermodynamics and kinetics of ethylene carbonate reduction at the electrode–electrolyte interface with molecular simulations. J. Chem. Phys. 155, 204703 (2021).

58. Appel, A. M. & Helm, M. L. Determining the overpotential for a molecular electrocatalyst. ACS Catal. 4, 630–633 (2014).

    Article  CAS  Google Scholar 

59. Rudshteyn, B. et al. Calculation of metallocene ionization potentials via auxiliary field quantum Monte Carlo: toward benchmark quantum chemistry for transition metals. J. Chem. Theory Comput. 18, 2845–2862 (2022).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

60. Laoire, C. O., Plichta, E., Hendrickson, M., Mukerjee, S. & Abraham, K. M. Electrochemical studies of ferrocene in a lithium ion conducting organic carbonate electrolyte. Electrochim. Acta 54, 6560–6564 (2009).

    Article  CAS  Google Scholar 

61. Hall, D. S., Self, J. & Dahn, J. R. Dielectric constants for quantum chemistry and Li-ion batteries: solvent blends of ethylene carbonate and ethyl methyl carbonate. J. Phys. Chem. C 119, 22322–22330 (2015).

    Article  CAS  Google Scholar 

62. Spotte-Smith, E. W. C. et al. Quantum chemical calculations of lithium-ion battery electrolyte and interphase species. Sci. Data 8, 203 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

63. Wang, Y., Nakamura, S., Ue, M. & Balbuena, P. B. Theoretical studies to understand surface chemistry on carbon anodes for lithium-ion batteries: reduction mechanisms of ethylene carbonate. J. Am. Chem. Soc. 123, 11708–11718 (2001).

    Article  CAS  PubMed  ADS  Google Scholar 

Download references

Author notes

1. These authors contributed equally: James M. Stevenson, Garvit Agarwal, Sibali Debnath, Zachary K. Goldsmith, Monosij Mondal, Asya Svirinovsky-Arbeli, Richard A. Friesner

Authors and Affiliations

1. Schrödinger, New York, NY, USA

   James M. Stevenson, Garvit Agarwal, Rishabh D. Guha & John L. Weber

2. Department of Chemistry, Columbia University, New York, NY, USA

   Sibali Debnath, Zachary K. Goldsmith, Monosij Mondal & Richard A. Friesner

3. Department of Chemical Engineering, Columbia University, New York, NY, USA

   Asya Svirinovsky-Arbeli, In Won Yeu, Amy S. Metlay, Yongbeom Kwon, Lakshmi Bhai, Daniel A. Steingart, Lauren E. Marbella & Alexander Urban

4. New York Structural Biology Center, New York, NY, USA

   Eric G. Keeler

5. Schrödinger, Portland, OR, USA

   Leif D. Jacobson

6. Schrödinger, San Diego, CA, USA

   Mathew D. Halls

Contributions

J.M.S. conceptualized the project, developed the MLFF methodology, curated data, investigated predictions, validated the model and authored the paper. G.A. curated data, developed the MLFF methodology, performed molecular dynamics simulations, analysed the results and edited the paper. S.D. conceptualized and performed NMR predictions, analysed the results and edited the paper. Z.K.G developed clustering methodology, validated the MLFF model, developed methodology for and performed reduction potential calculations, analysed the results, and authored and edited the paper. M.M. developed the clustering methodology, conceptualized and performed NMR simulations, validated the MLFF model, analysed the results, and authored and edited the paper. A.S.-A. developed the NMR methodology, collected and analysed the NMR spectra, and wrote and edited the paper. I.W.Y. performed AIMD simulations, analysed results and edited the paper. A.S.M. performed solution NMR experiments, density and viscosity measurements, analysed results and edited the paper. Y.K. performed Raman measurements, analysed results and edited the paper. E.G.K. helped set up the NMR measurements and edited the paper. L.B. performed the NMR experiments, analysed results and wrote the paper R.D.G. performed validation of our MLFF and other force fields and edited the paper. J.L.W. analysed the results and edited the paper. D.A.S. authored and edited the paper. L.E.M. supervised the project, developed the NMR methodology, analysed the results, and authored and edited the paper. A.U. supervised the project and edited the paper. L.D.J. developed the MLFF methodology, conceptualized the project and edited the paper. M.D.H. supervised the project. R.A.F. conceptualized the project, supervised the project, authored and edited the paper, and acquired funding.

Corresponding authors

Correspondence to John L. Weber, Lauren E. Marbella or Richard A. Friesner.

Competing interests

R.A.F. has a considerable financial stake in, is a consultant for, and is on the scientific advisory board of Schrödinger. J.M.S. and L.D.J. were employed by Schrödinger at the time of this writing. G.A., R.D.G., J.L.W. and M.D.H. are employed by Schrödinger. The other authors declare no competing interests.

Peer review information

Nature thanks Steven Greenbaum, who co-reviewed with Giselle de Araujo Lima e Souza; and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions