Excess entropy scaling explains the enhanced dynamics of the ionic liquid 1-ethyl-3-methylimidazolium chloride in external electric fields

Abstract

In this study, we test scaling relationships between excess entropy and self-diffusion of the ionic liquid (IL) 1-ethyl-3-methylimidazolium chloride ([C₂C₁im][Cl]) across temperatures and varying external electric fields (EEFs). EEFs are observed to increase the dynamics of the IL by altering the structure of the liquid. This study assesses the validity of excess entropy scaling relationships for predicting self-diffusion coefficients of an IL across various temperatures and EEFs using molecular dynamics simulations. We compute the excess entropy of [C₂C₁im][Cl] to quantify EEF-induced structural changes and diffusion behavior. Results show that EEFs primarily affect anion–anion organization at short and medium length scales, while cation–cation ordering shifts only at long length scales. Additionally, counterion configurational ordering dominates excess entropy contributions over radial or translational organization. These findings establish excess entropy scaling as a robust tool for describing IL transport under EEFs, informing the design of solvents for sustainable separations under external stimuli.

---

Introduction

Many modern technologies, such as permanent magnets, batteries, and catalysts rely on rare earth elements (REEs). Chemical separation processes are key to producing and recycling REEs.¹ A 2019 National Academies report highlighted that understanding of external stimuli that result in physical changes could lead to advancements in chemical separations.² Ionic liquids (ILs) offer the opportunity to develop this understanding given that they are composed of ions which respond to an external electric field (EEF). ILs are a class of solvent extensively studied for separation processes; they are liquid below 100 °C, and offer high thermal stability, low volatility, and tunability of their thermophysical properties.³ Recent studies have explored the influence of EEFs on the thermophysical properties of ILs, finding direct coupled effects on dynamic and structural properties.^4–7 These EEF-induced property changes could be exploited to improve the performance of ILs in REE separation processes. However, relationships used to describe the coupled effect on the structure and dynamics of ILs across EEFs and temperature ranges have not been tested. Here we use a combination of molecular dynamics (MD) simulation and excess entropy scaling relationships to describe self-diffusion coefficients of ILs across various temperatures and EEFs.

Excess entropy (S_exc) scaling relationships have been used to describe the transport properties of multiple chemical systems across many thermodynamic states.^8–12 For a given density and temperature, S_exc(ρ, T) quantifies the deviation of a system's entropy S(ρ, T) from the entropy of an ideal gas

S_exc(ρ, T) ≡ S(ρ, T) − S_id(ρ, T). (1)

Rosenfeld introduced S_exc scaling relationships first by showing the entropy dependence of self-diffusion coefficients for hard-sphere, soft-sphere, Lennard-Jones, and single-component plasma systems with MD simulations.¹³ Later, Dzugutov elucidated this relationship with MD simulations, by demonstrating a direct connection between structure and diffusion in a dense liquid.¹⁴ Dzugutov argued that the number of accessible configurations, quantified by S_exc, is proportional to the frequency of local structural relaxations, which defines the rate of diffusion across a medium.¹⁴ This allowed Dzugutov to establish a scaling relationship between S_exc and diffusion coefficients, and demonstrated the relationship by showing how reduced diffusion coefficients can be predicted by estimating the number of accessible configurations from structural correlations obtained from MD simulations which S_exc(ρ, T) quantifies. For an expanded theoretical and analytical overview of excess entropy scaling, we refer the reader to Dyre's comprehensive perspective on the subject.⁸ Excess entropy scaling relationships have also been explored for systems like water, beryllium fluoride, germanium oxide, silicon dioxide, silicon, and germanium.^15–18 Lötgering-Lin and Gross proposed an entropy-scaling approach for Newtonian shear viscosities based on group contributions and demonstrated examples of S_exc scaling relationships for n-alkanes and ketones.¹⁹ Wang and Ghaffarizadeh also demonstrated that this kind of relationship applies to complex fluids like active matter systems, and leveraged it to understand how the transport properties of an inactive fluid react to the addition of active particles that can convert free energy into motion.²⁰

S _exc scaling relationships have also been tested for ILs such as dimethylimidazolium chloride. Malvaldi employed MD simulations to investigate if the Rosenfeld and Dzugutov S_exc scaling laws could describe the diffusion coefficients of ILs.²¹ They observed large deviations from Rosenfeld scaling when only translational degrees of freedom were used to estimate S_exc. Only after computing S_exc by considering orientational and configurational degrees of freedom, in addition to translational, did the diffusion coefficients display Rosenfeld scaling. Malvaldi also observed that diffusion coefficients versus S_exc for the anion and the cation superimposed up to the deeply supercooled regime. In addition, they presented evidence for Dzugutov scaling laws being able to predict the ratio between cation and anion mobility by capturing the higher mobility of the heavier and larger cation typically observed in imidazolium-based ILs.²² This is of particular interest since it shows that S_exc scaling relationships can bridge structure and dynamics for systems with complex dynamics like ILs where different ions experience different scales of mobilities and nanostructure, just as other complex systems like binary glasses and active matter systems do.^23,24

Motivated by the observed increases in diffusion coefficients of ILs in the presence of EEFs and the need to describe transport coefficients in external stimuli, we test for the generalizability of excess entropy scaling relationships for the IL 1-ethyl-3-methylimidazolium chloride [C₂C₁im][Cl] in EEFs.^4,25 We use MD simulations of [C₂C₁im][Cl] to compute excess entropies across EEFs and temperature ranges and test if computed self-diffusion coefficients from the simulation scale with excess entropy computed from the IL structure.

Methods

Simulation details and protocol

The large-scale atomic/molecular massively parallel simulator (LAMMPS) package was used for all MD simulations along with the general Amber force field (GAFF) to describe inter- and intramolecular interactions.^26–28 Partial atomic charges were computed using the Gaussian 09 package and the restrained electrostatic potential (RESP) method, based on the optimized vacuum structure of each ion at the B3LYP/6-311++G(d,p) level.²⁹ Then, all charges were uniformly scaled by 0.8 to account in an approximate way for polarization effects.^30,31 The particle–particle particle–mesh solver (PPPM) was used to compute long-range Coulombic interactions with a 12 Å real space cutoff. A cutoff of 12 Å was used for Lennard-Jones interactions along with long-range tail corrections to the energy and pressure. For initial configurations, 1007 pairs of [C₂C₁im][Cl] were randomly placed in a cubic box using PackMol.³² This was followed by an energy minimization for 10 000 iterations or until one of the stopping criteria was met (a tolerance of 1 × 10⁻⁴ kcal mol⁻¹ for energy and 1 × 10⁻⁶ kcal mol⁻¹ Å⁻¹ for force). The temperature damping parameter was 100 fs for both canonical (NVT) and isothermal–isobaric (NPT) ensembles, while the pressure damping parameter for the NPT ensemble was 1000 fs. A time step of 1 fs was used for all simulations. All simulations were carried out at a pressure of 1 atm. LAMMPS input files with detailed specifications for pair styles and force field parameters can be found in the SI and Zenodo repository.

Unfortunately, to the best of our knowledge, there are no experimental measurements of self-diffusion coefficients for any ionic liquid under external electric fields. Thus, our simulations are purely predictive at this point. However, we have validated the force field used for the no external electric field case. The computed density for [C₂C₁im][Cl] at 360 K with no EEF was 1108.5 kg m⁻³ which deviates less than 1% from the closest experimental measurement at 359 K by Součková et al. of 1109.1 ± 1.1 kg m⁻³.³³ Therefore, we believe the force field describes the structure of the IL since it captures the experimental packing well. In addition, another computational study employing another force field is available to compare the diffusion coefficient for both ions. For the present study at 360 K with no EEF, the computed self-diffusion coefficient for the [C₂C₁im] ion was 6.73 × 10⁻¹¹ m² s⁻¹ while for [Cl] it was 6.64 × 10⁻¹¹ m² s⁻¹. This compares well to the values computed by Jiang et al. at 353 K, (5.925 ± 10.06) × 10⁻¹¹ m² s⁻¹ and (4.365 ± 8.11) × 10⁻¹¹ m² s⁻¹ respectively.³⁴

While polarizable models may be more accurate at describing the effect of temperature on thermodynamic and dynamic properties, we find that the force field used for this study also captures these trends well enough to demonstrate the relationship between diffusion and excess entropy. Specifically, this force field is capable of capturing trends across diffusion and viscosity for a range of ILs, which is challenging to capture.^35,36 In addition, the force field captures the structure of the IL which is critical for estimating excess entropy. Thus, we anticipate these trends would translate well enough to test excess entropy scaling relationships for this system. A future direction may be to explore if a polarizable force field renders a stronger excess entropy relationship than a non-polarizable force field.^37,38

The simulation protocol was the following. After minimizing initial configurations, a static and homogeneous EEF of varying strengths (0.0000 to 1 × 10⁻¹ V Å⁻¹) was applied only along the direction of the x-axis for all simulations. The lower bound of EEF selected in this study was defined from other studies where changes in structure are observed only until EEFs of 2 × 10⁻² V Å⁻¹ or larger are applied.^39,40 It is important to note that the magnitudes of EEF used are necessary to observe effects within the time scales of MD simulations, however they may not be experimentally achievable. The upper limit is defined from the stability of the thermostats where for EEF larger than 1 × 10⁻¹ V Å⁻¹ artifacts may be introduced by the thermostat.⁴¹ Careful application of an MD thermostat in EEF is necessary to avoid introducing artifacts through the thermostat that lead to the violation of equipartition.⁴² Lan and Wang characterized such artifacts and quantified the violation of equipartition of multiple thermostats in EEFs by computing the parameter .⁴³T_x and T_y are ensemble average temperatures computed from the velocity of atoms along the x and y directions respectively. T_x corresponds to the temperature in the direction of the EEF. When χ approaches 1, violation of equipartition is strong meanwhile when χ approaches 0 there is little to no violation of equipartition. We have previously shown that the magnitudes of EEF and the methods used to study similar ILs are appropriate by computing small χ values that approach 0 rather than 1.⁷ After an EEF was applied, an annealing NPT step at 750 K for 1 ns was executed. Then, the system was cooled to the target temperature over 4 ns. After reaching the target temperature, the equilibrium ensemble density was computed from a 7 ns trajectory outputting box density every 1000 fs. After an equilibrium density ensemble average was obtained, the box dimensions were uniformly deformed according to the obtained equilibrium density, and a 5 ns dynamics and structure equilibration step was performed in the NVT ensemble. Finally, a 20 ns NVT dynamics production run was executed where the positions of all atoms were saved every 10 ps to compute the self-diffusion coefficient of each ion, accounting for the drift induced by the EEF.²⁵ These positions were also used to compute radial, configurational, and spatial distribution functions.

Radial, configurational and spatial pair distribution functions

From the stored positions of the NVT production run, radial distribution functions g_ij(r) between the center of mass of ions i and j were computed with the package PyLAT using a 0.1 Å resolution.⁴⁴ Configurational pair distribution functions g^ij_conf(ω₁|r), where ω₁ = {−π < ϕ < π, 0 < θ < π} describes the relative spatial orientation around ions, were constructed from spatial distribution functions obtained through the TRAVIS trajectory analyzer.⁴⁵ We computed a spatial distribution up to a maximum observation radius of 28 Å with a bin number in all directions of 300 in reference to the cation and with the reference axis shown in Fig. 1. Then, the polar coordinates of each voxel in space were computed in reference to the axis shown in Fig. 1.

Fig. 1 (Left) Cation reference axis is depicted with red, green, and blue arrows representing the positive x, y, and z axes. Colored rings of radius 4.5 Å serve as visual aids: the red ring indicates variations in ϕ, and the blue and green rings illustrate variations in θ. (Right) Spherical reference axis showing how θ and ϕ vary. Intersections of the colored rings correspond to specific θ and ϕ values—for example, where green and blue rings meet in the positive z-direction θ = 0 and where red and blue rings meet in the positive x-direction θ = π/2 and ϕ = 0. The point where green and red rings meet in the positive y-direction corresponds to θ = π/2 and ϕ = π/2, while the point where the rings meet again in the negative y direction correspond to θ = π/2 and ϕ = −π/2.

Once the polar coordinates r, ϕ and θ were obtained, histograms were constructed and values of spatial relative density obtained from TRAVIS were projected across ϕ and θ for a given r from 0 to 28 Å with a resolution of 0.1 Å.

Self-diffusion coefficients

We computed the self-diffusion coefficients of each ion in the presence of EEFs with the method presented in our previous work.²⁵ We computed the mean square displacement (MSD) of the center of mass of each ion for each time step with the saved positions from the NVT production run. Using the MSD of each ion, the Einstein relationship was used to compute the self-diffusion coefficient.⁴⁶ These computations were done with a modified branch of PyLAT that considers the drift induced by the EEF.^25,44 The uncertainty in diffusion coefficients was estimated by computing the standard deviation between three independent replicates of [C₂C₁im][Cl] with different starting positions and velocities.

Scaling relationships and reduced units

In general, the scaling relationship between S_exc and the diffusion coefficient takes the formwhere c₁ and c₂ are material-specific constants, D* is the reduced diffusion coefficient, and is the per-particle excess entropy for a system with N particles.²⁰ For this study N is the total number of ions defined as N = N_cat + N_ani where N_cat is the number of cations and N_ani is the number of anions. For pure systems, D* is reduced with macroscopic thermodynamic observables like temperature and density or microscopic quantities like collision frequencies and diameters. For ILs which are a mixture of cations and anions, D* is reduced as a combined diffusion coefficient in the following wayIn eqn (3), x is the fraction of the total ions defined as x = N_cat/N or x = (N − N_ani)/N. D_i={cat,ani} is the self-diffusion coefficient of the anion or the cation. This corresponds to Dzugutov scaling where χ_i={cat,ani} is a microscopic scaling factor that considers the same species' collision frequencies and diameters along with their cross-term counterparts. χ_i={cat,ani} is computed in the following way:where σ_ij is the kinematic radius, which is defined as the distance r where the highest peak of the radial distribution function max(g_ij(r)) is observed. ρ_j is the number density of the cation or anion. m_i and m_j correspond to the molar masses of the ions of the IL.

Excess entropy estimations

Excess entropy estimates all the possible configurations available to a system at a given density and temperature. Derivations show that this quantity can be estimated through pair correlation functions.⁸ For single-component systems with spherical molecules, this can be completely captured using only translational pair distribution functions, better known as radial distribution functions. For mixtures of non-spherical molecules, translational pair distribution functions average over the configurational ordering that ions have around each other and do not include orientational degrees of freedom. To completely estimate excess entropy for an IL, translational, configurational, and orientational degrees of freedom must be accounted for by computing translational, configurational, and orientational pair distribution functions. We will now describe how the excess entropy of the IL was estimated for each temperature and EEF.

The excess entropy for the IL is expressed as a partial molar excess entropy, computed from the excess entropy in reference to the cation and the anion with the following equation

In eqn (5), x is related to the IL stoichiometry, which in this work is 0.5. The excess entropy in reference to each ion is computed with the following equationwhere g^ij₍₂₎(r, ω₁, ω₂) is the two-particle correlation function in terms of the distance r between the center of mass of ions i and j, and the spatial orientation around each ion and their relative alignment.⁴⁷ The relative spatial orientation around ions in g^ij₍₂₎(r, ω₁, ω₂) is described by the polar angle pair ω₁ = {ϕ, θ} with a reference axis defined by the reference ion i in eqn (6). The range for ϕ is from −π to π and for θ from 0 to π. The relative alignment between ions is described by the Euler angles ω₂ = {α, β, γ}. The ranges for α, β, and γ are from 0 to π, 0 to 2π, and 0 to 2π, respectively. ρ is the system's number density and Ω = 8π² is a normalization constant. In eqn (6)ρ is the number density. Therefore, the indefinite integral in eqn (6) is resolved as,To our benefit, the two-particle correlation function can be factored into its translational, configurational, and orientational components as

g^i,j₍₂₎(r, ω₁, ω₂) = g^i,j_trans(r)·g^i,j_conf(ω₁|r)·g^i,j_orient(ω₂|ω₁, r) (8)

where g^i,j_trans(r) is the typical g_i,j(r) and describes translational correlations.⁴⁸ The relative spatial orientation around each ion is described by the configurational distribution function g^i,j_conf(ω₁|r). g^i,j_conf(ω₁|r) is an extension of g^i,j_trans(r) in which for a given radius r, g^i,j_conf(ω₁|r) is a surface describing the relative bulk density of ion j with respect to a reference ion i across polar angles ω₁. g^i,j_orient(ω₂|ω₁, r) describes the orientational component, which for a given r and combination of ω₁ angles describe the relative bulk density of ions j orienting themselves with a reference ion i across Euler angles ω₂. The factorization of g^i,j₍₂₎(r, ω₁, ω₂) allows us to combine eqn (8) and (6) and distribute the integrals to factor out each translational, configurational and orientational contribution to the excess entropy as

Mavaldi showed that for dimethylimidazolium chloride ([C₁C₁im][Cl]), the translational and configurational contributions can provide a reasonable estimate of total excess entropy to obtain a Dzugutov type scaling relationship.²¹ Because [Cl] is spherical and does not have a molecular frame of reference, we don’t consider configurational or orientational contributions in reference to it and only compute translational contributions as

The translational contributions are computed asFollowing the same approach as Mavaldi, we estimate the excess entropy contribution from the cation [C₂C₁im] considering only translational and configurational contributions asThe configurational contributions are computed as

We validated our approach for computing S_exc by comparing our estimations for 3 replicates for a SPC water system with the results that Zielkiewicz obtained when measuring the orientational contributions to the S_exc of water.⁴⁷ We find that our methods estimate the configurational excess entropy well; we estimate around −17.70 ± 0.06 J mol⁻¹ K⁻¹. This compares well with the value of −14 J mol⁻¹ K⁻¹ computed by Zielkiewicz. Further details and information can be found in the SI.

Results and discussion

Structure and dynamics are coupled in an EEF

The counterion center of mass spatial distribution function in Fig. 2 shows the effect of the EEF on the structure of [C₂C₁im][Cl] for replicate 1 as a representative example. For the zero EEF case, the [Cl] ion localizes around the imidazolium ring, consistent with the charge alternating structure present in ILs.⁴⁹ There is a depletion zone around the non-polar ethyl group of the cation. For the 1 × 10⁻¹ V Å⁻¹ EEF case, the [Cl] ions populate similar positions as for the zero EEF case, but the iso-surface without an EEF is larger than the iso-surface with an EEF of 1 × 10⁻¹ V Å⁻¹. This indicates that the EEF disrupts [C₂C₁im] and [Cl] spatial correlations, causing [Cl] to localize around [C₂C₁im] more than when there is no EEF. Fig. 3 shows the self-diffusion coefficients for each ion at multiple temperatures and across multiple EEFs. For all temperatures, the presence of a strong enough EEF increases the self-diffusion coefficient for both ions.

Fig. 2 Spatial distribution function of the [Cl] anion around the center of mass of the [C₂C₁im] cation showing an iso-surface with a cutoff 3 times the bulk density. The left panel shows the spatial distribution function without an EEF. The right panel shows the spatial distribution when an EEF of 1 × 10⁻¹ V Å⁻¹ is present.

Fig. 3 Self-diffusion coefficients for [C₂C₁im] and [Cl] on the left and right panel respectively for multiple temperatures across different EEFs. At high EEF there is a substantial increase in the self-diffusion coefficient for both ions. Both panels share the same legend. For both ions, self-diffusion coefficients are unaffected and overlap for EEF up to 2 × 10⁻² V Å⁻¹. Only when an EEF of 4 × 10⁻² V Å⁻¹ or larger is applied are changes in self-diffusion coefficients observed. Error bars represent the standard deviation between three independent replicates of [C₂C₁im][Cl].

Fig. 2 and 3 suggest that the structural changes caused by the EEF are coupled with the change in dynamics of [C₂C₁im][Cl]. Therefore, we hypothesize that the effect on dynamics could be understood through the effect on liquid structure or vice versa. We will now show that excess entropy scaling relationships can model such behavior, allowing us to map IL dynamics to liquid structure across multiple thermodynamic states, including EEFs.

Translational excess entropy describes the effect of EEF on short and medium range IL structures

The translational excess entropy from the [C₂C₁im]–[C₂C₁im], [C₂C₁im]–[Cl] and [Cl]–[Cl] translational correlations was computed across multiple simulations as a function of EEF strength to quantify the effect of EEFs on IL structure and to test excess entropy scaling relationships. The translational excess entropy contribution to the total excess entropy from the g(r) was computed with eqn (11). Panels A, B, and C of Fig. 4 show the g(r) between [C₂C₁im]–[C₂C₁im], [C₂C₁im]–[Cl] and [Cl]–[Cl] at 380 K across EEFs for replicate 1 as a representative example. Like self-diffusion coefficients, once a high enough EEF is applied, an effect on IL structure is observed. For the g(r) between [C₂C₁im] and [C₂C₁im] in Fig. 4 panel A, a shoulder around 4.5 Å was observed to disappear at 1 × 10⁻¹ V Å⁻¹. In addition, the g(r) between [C₂C₁im] and [C₂C₁im] at 1 × 10⁻¹ V Å⁻¹ shows less structure across r compared to the g(r) at lower EEFs. Specifically, the maximum of g(r) decreases and the minimum of g(r) around 12 Å increases. This is interpreted as a loss in correlation in the [C₂C₁im] and [C₂C₁im] structure, which would correspond to a less negative excess entropy, corresponding to an increase in mobility. A similar effect is observed for the [C₂C₁im]–[Cl] and [Cl]–[Cl] g(r) where the maximum of g(r) decreases and minimum of g(r) in valleys increases. Specifically, for the g(r) between [Cl] and [Cl], peaks and valleys are smoother and wider with local maxima and minima also shifting to longer distances. These observations are direct effects of the EEF on IL structure, with ion and counterion radial organization being perturbed by the EEF, which affects the liquid packing. For example, the first peak in [Cl]–[Cl] g(r) decreases and is shifted to a longer distance, suggesting a larger volume. These observations are consistent with molar volume increasing with an increasing EEF.⁷ Fig. S5–S8 in the SI demonstrate the effect of the EEF to features of the g(r) in clearer detail and compares them to g(r) processed with KAMEL-LOBE, a technique used to reduce noise in g(r).⁵⁰

Fig. 4 Radial distribution functions between [C₂C₁im]–[C₂C₁im], [C₂C₁im]–[Cl] and [Cl]–[Cl] on panels A–C, respectively, as a function of EEF at 380 K. Estimates of translation excess entropy in [J mol⁻¹ K⁻¹] as a function of radius from [C₂C₁im]–[C₂C₁im], [C₂C₁im]–[Cl], and [Cl]–[Cl] structure on panels D–F, respectively, for multiple EEFs at 380 K computed with eqn (11).

Translational excess entropy can be used to describe how the EEF affects the liquid structure of the IL across r. The translational excess entropy from the [C₂C₁im]–[C₂C₁im], [C₂C₁im]–[Cl] and [Cl]–[Cl] translational correlations computed from each respective g(r) using eqn (11) are shown in Fig. 4 panels D, E, and F. In panel D of Fig. 4 it can be observed that the EEF affects the medium range structure of [C₂C₁im]–[C₂C₁im] more than the short-range structure since the translational excess entropy across r starts to deviate after 10 Å. Beyond 10 Å, as the magnitude of the EEF increases, the translational excess entropy from [C₂C₁im]–[C₂C₁im] correlations become less negative. This means that the EEF effects [C₂C₁im]–[C₂C₁im] structure by reducing medium range structural correlations. Panel E of Fig. 4 shows the translational excess entropy from the [C₂C₁im]–[Cl] translational correlations. The EEF affects the short and mid-range translational correlations between [C₂C₁im] and [Cl] starting at 4.5 Å. Then, the mid-range [C₂C₁im]–[Cl] translational correlations between 5 Å and 10 Å are affected. Finally, large effects on translational correlations from the EEF are observed after 10 Å. Panel F of Fig. 4 also shows that the EEF affects the short and medium range of the [Cl]–[Cl] translational correlations starting at 4.5 Å. Beyond 4.5 Å, the translational excess entropy from the [Cl]–[Cl] translational correlations are less negative as the EEF increases. But most importantly, the EEF alters the [Cl]–[Cl] translational excess entropy more than the [C₂C₁im]–[C₂C₁im] translational excess entropy. The changes in structural excess entropy for [Cl]–[Cl] compared to [C₂C₁im]–[C₂C₁im] suggest that the anion has a larger role than the cation in the effect that an EEF has on changing thermophysical properties of the IL. The role of the anion in response to the electric field may be due to its higher charge density compared to the cation, given that the charge density of ions is correlated to the response of the IL to the EEF.⁷ However, the largest change in translational excess entropy is observed between [C₂C₁im]–[Cl] correlations, which Fig. 2 also demonstrates by the differences in the size of the regions where the [Cl] would localize. The decrease in size of the region where the [Cl] is localized suggests more disorder and a less negative excess entropy which correspond to higher dynamics.

Translational excess entropy alone cannot describe the differences in self-diffusion coefficients between [C₂C₁im] and [Cl]. By definition S_exc(ρ, T) ≡ S(ρ, T) − S_id(ρ, T) quantifies how different the real entropy of the system is compared to an ideal gas at the same temperature and density. S_exc(ρ, T) will always be negative because for a given density and temperature, S_id(ρ, T) would always be larger than S(ρ, T). The more negative S_exc(ρ, T) is, the smaller S(ρ, T) would be compared to S_id(ρ, T). Per Boltzmann's equation for entropy S(ρ, T) = k_B ln Ω(ρ, T), smaller values of S(ρ, T) correspond to a smaller Ω(ρ, T) which quantifies the available microstates for a given macrostate. Excess entropy scaling relationships hypothesize that there is an indirect proportionality between the transport properties of a system and Ω(ρ, T). That is, a system with a small Ω(ρ, T) would have smaller self-diffusion compared to that of a system with a larger Ω(ρ, T). In terms of S_exc(ρ, T), excess entropy scaling relationships predict that the more negative S_exc(ρ, T) is, the smaller the self-diffusion coefficient would be. The computed translational excess entropy from the [C₂C₁im]–[C₂C₁im] and [Cl]–[Cl] translational correlations alone do not follow this proportionality. From panels D and F of Fig. 4, the translational excess entropy from [C₂C₁im]–[C₂C₁im] translational correlations converged to less negative values than the translational excess entropy from [Cl]–[Cl]. According to excess entropy scaling relationships, a less negative translational excess entropy would correspond to higher self-diffusion coefficients for [C₂C₁im] compared to that of [Cl] which is not observed in Fig. 3. Malvadi also observed that translational excess entropy alone cannot be used for excess entropy scaling relationships for the diffusion coefficients of dimethyl imidazolium chloride.²¹

Configurational excess entropy describes and quantifies the effect of EEF on long range IL structure

The configurational distribution function g_conf(r|ω₁) between [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] was computed across EEFs and temperatures to understand and quantify the effect of EEF on IL structure. g_conf(r|ω₁) with the [Cl] ion as a reference was not computed because [Cl] does not have a molecular frame of reference. In g_conf(r|ω₁) the pair ω₁ = {ϕ, θ} denotes the spherical coordinates with distance r of a point in space with the reference axis shown at the center of Fig. 5. g_conf(r|ω₁) was computed by following the structural analysis protocol provided in the Methods section and was used to estimate the configurational excess entropy with eqn (13). g_conf(r|ω₁) is an extension of g(r) and holds information on how ions organize around each other at a given r. In a typical g(r), the spatial distribution of species is averaged and projected onto one dimension r. Therefore, a typical g(r) can only describe variations in spatial density along r. In contrast, g_conf(r|ω₁) can describe spatial variations for a fixed r across spherical coordinates ω₁ with ω₁ = {−π < ϕ < π, 0 < θ < π}. Therefore, g_conf(r|ω₁) is a surface that describes how the spatial relative density of [C₂C₁im] center of mass (COM) or [Cl] varies at a fixed r around the COM of a reference [C₂C₁im] ion.

Fig. 5 Configurational distribution function g_conf(r|ω₁) between [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] at 380 K with zero EEF in panels A and B, respectively, for r = 4.5 Å. Panels C and D show g_conf(r|ω₁) between [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] at 380 K with an EEF of 1 × 10⁻¹ V Å⁻¹ at the same radius. The top right color bar is the heatmap for panels A and C while the bottom color map is the heatmap for panels B and D. In g_conf(r|ω₁), r and ω₁ = {−π < ϕ < π, 0 < θ < π} are spherical coordinates within the reference axis shown in the center panel. The origin is the [C₂C₁im] COM and r is the distance from the origin to another [C₂C₁im] COM or [Cl]. The radius r = 4.5 Å corresponds to the maximum peak in the [C₂C₁im]–[Cl] g(r) and the shoulder in the [C₂C₁im]–[C₂C₁im] g(r). In the middle the cation reference axis is shown with red, green, and blue arrows representing the positive x, y, and z axes. Colored rings of radius 4.5 Å serve as visual aids: the red ring indicates variations in ϕ, and the blue and green rings illustrate variations in θ. The reference axis in the middle right panel are visual aids showing how θ and ϕ vary. The reference axis defined and described in Fig. 1.

Each panel in Fig. 5 describes how the COM of [C₂C₁im] or [Cl] organizes around a reference [C₂C₁im] at a temperature of 380 K with and without an EEF for replicate 1 as a representative example. The reference [C₂C₁im] molecule is shown in the center panel of Fig. 5. The reference axis origin is at the [C₂C₁im] COM. A description of the reference axis is provided in the caption of Fig. 1. Rings in the center panel serve as visual guides and are placed at the origin with a radius of 4.5 Å. g_conf(r|ω₁) is specifically shown for r = 4.5 Å since it corresponds to the position of the maximum peak and shoulder in the [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] g(r), respectively. Panel A in Fig. 5 shows g_conf(r = 4.5 Å|ω₁) between [C₂C₁im]–[Cl] without an EEF. In panel A only peaks were observed at and along various angles for ϕ, indicating that [Cl] aggregates mostly in the plane of the imidazolium ring. That is, [Cl] ions mostly aggregate along the red ring around the reference [C₂C₁im] ion in the center panel of Fig. 5. Specifically, [Cl] ions aggregate in three directions along the red ring, since 3 distinct peaks across ϕ and along are observed. The wide medium peak around ϕ = 0 and corresponds to [Cl] aggregating near the direction of the C2 atom in between the two N atoms in the imidazolium ring. As ϕ gets closer to π along , a large sharp peak is observed. This peak corresponds to [Cl] aggregating near the C5 atom. As values for ϕ go from 0 to −π, a valley is observed followed by a small peak in the direction of the C4 atom. Finally, at ϕ = −π the continuation of the peak in the direction of ϕ = π and was observed. In summary, the [Cl] anions align along the plane of the imidazolium ring and aggregate in the directions of the C atoms in the imidazolium ring, consistent with the results shown in Fig. 2 and 5.

Fig. 5 panel B shows the g_conf(r = 4.5 Å|ω₁) between [C₂C₁im]–[C₂C₁im] without an EEF. In panel B only the peaks near θ = 0 and θ = π are observed along various angles for ϕ. This indicates that [C₂C₁im] aggregates mostly in the directions where the blue and green rings meet, meaning that at this radius, other [C₂C₁im] ions aggregate normal to the imidazolium ring and mostly avoid the red ring shown in the center of the panel. However, there are small peaks along that show that there are directions along the red ring that [C₂C₁im] aggregate. Specifically, this occurs for ϕ near zero and −π, which corresponds to the direction where the red and blue ring meet in the positive and negative x directions, respectively. These points correspond to the direction of the C2 atom in the imidazolium ring between C4 and C5. These observations show how g_conf(r = 4.5 Å|ω₁) can provide further detail to the nature of the IL structure and allow us to understand the effects of EEF on IL nanostructure that result in changes to their thermophysical properties.

The effect of EEFs on [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] configurations are observed by comparing panel A with panel C and panel B with panel D in Fig. 5. At r = 4.5 Å no drastic changes in configurations are observed since there are no shifts in the positions of the peaks across any g_conf(r = 4.5 Å|ω₁) surface. However, there is an overall decrease in relative density in g_conf(r|ω₁) for [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im]. Panel C shows g_conf(r = 4.5 Å|ω₁) for [C₂C₁im]–[Cl] in the presence of an EEF. The position of the peaks does not change, however the height of each peak does. Specifically, the height of the peaks at ϕ = 0 and ϕ = π are lowered. The peak at ϕ = 0 is narrower across θ compared to the same peak for g_conf(r = 4.5 Å|ω₁) computed without an EEF. In panel D, the same effect was observed, where the heights of the peaks for g_conf(r = 4.5 Å|ω₁) between [C₂C₁im]–[C₂C₁im] decreased. That is, for both [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] configurations, the EEF seems to smooth peaks and lower relative density values for g_conf(r = 4.5 Å|ω₁). Thus, the EEF decreases spatial correlations between counterions.

As previously mentioned, we computed g_conf(r|ω₁) to quantify the effect of EEFs on IL liquid structure and test excess entropy scaling relationships for the self-diffusion coefficients of ILs in EEFs. To test excess entropy scaling relationships, we compute the configurational excess entropy from the [C₂C₁im]–[Cl] g_conf(r|ω₁) and [C₂C₁im]–[C₂C₁im] g_conf(r|ω₁) with eqn (13). In Fig. 6, we show the configurational excess entropy as a function of r for replicate 1 as a representative example. We observe a clear effect of the EEF on the configurational excess entropy, especially on the [C₂C₁im]–[Cl] configurational excess entropy. Like the translational excess entropy in Fig. 4, we observe that the configurational excess entropy is more negative as r increases and less negative as the EEF increases.

Fig. 6 Estimate of configurational excess entropy in [J mol⁻¹ K⁻¹] as a function of radius from the [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] structure at 380 K across EEFs for replicate 1 as a representative example. The solid and dot dashed lines represent the configurational excess entropy between [C₂C₁im]–[C₂C₁im] and [C₂C₁im]–[Cl] computed from g_conf(r|ω₁) respectively with eqn (13). The configurational excess entropy contribution from the [C₂C₁im]–[C₂C₁im] configurational correlations is smaller and less sensitive to the EEF compared to the configurational excess entropy from the [C₂C₁im]–[Cl] configurational correlations.

In contrast to the translational excess entropy, we observe that the configurational excess entropy oscillates at longer distances than the translational excess entropy. In Fig. 4, the translational excess entropy for [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] converges beyond 15 Å, whereas the configurational excess entropy for [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] still oscillates up to 28 Å. These oscillations reveal that the configurational ordering of the IL extends to a longer length scale than translational ordering. Such oscillations may be interpreted as a signal from the alternating counterion shells as distance from the reference ion increases. The oscillations also reveal that a bigger system may be needed to obtain a completely converged configurational excess entropy at a longer r.

In spite of configurational excess entropy oscillations, configurational excess entropy and g_conf(r|ω₁) describe and quantify the effect of EEF on IL liquid structure across distances and configurations. For values of r below 4.5 Å the configurational excess entropy is zero since there are no ions at those distances and no configurational correlations. As r increases, the configurational excess entropy between [C₂C₁im]–[Cl] dominates the configurational excess entropy from [C₂C₁im]–[C₂C₁im] starting at around r = 4.5 Å. Also, the short-range configurational correlations between [C₂C₁im]–[Cl] are more sensitive to the EEF than the [C₂C₁im]–[C₂C₁im] short range configurational correlations. Fig. 5 panels A and B describe why the configurational excess entropy from [C₂C₁im]–[Cl] is more negative that the configurational excess entropy from [C₂C₁im]–[C₂C₁im]. Fig. 5 panel A shows g_conf(r = 4.5 Å|ω₁) for [C₂C₁im]–[Cl] where g_conf(r|ω₁) ranges from 0 to 13. In comparison, Fig. 5 panel B shows g_conf(r = 4.5 Å|ω₁) for [C₂C₁im]–[C₂C₁im] where g_conf(r|ω₁) ranges from 0 to 2. The larger peak height from the g_conf(r = 4.5 Å|ω₁) for [C₂C₁im]–[Cl] corresponds to a higher correlation and therefore a larger contribution to the configurational excess entropy. The configurational excess entropy between [C₂C₁im]–[C₂C₁im] only starts to deviate across EEFs around r = 4.5 Å which corresponds to the maximum peak position in the [C₂C₁im]–[C₂C₁im] g(r). Finally, as the magnitude of the EEF increases, the configurational excess entropy for [C₂C₁im]–[Cl] and [C₂C₁im]–[C₂C₁im] becomes less negative. Specifically, the EEF impacts more the configurational excess entropy from [C₂C₁im]–[Cl] than the configurational excess entropy between [C₂C₁im]–[C₂C₁im].

Once the configurational excess entropy is accounted into the total excess entropy contribution from each ion, excess entropy can describe the differences between the diffusion coefficients of ions. Fig. S9 shows the total excess entropy from each ion computed with eqn (10) and (12) for [Cl] and [C₂C₁im] respectively for replicate 1. As the EEF is applied and the temperature increases, the total excess entropy from each ion is less negative. According to excess entropy scaling relationships, a less negative total excess entropy would correspond to higher self-diffusion coefficients for [Cl] compared to that of [C₂C₁im] which is observed in Fig. 3.

Excess entropy scaling test

We computed the partial molar excess entropy and the combined reduced self-diffusion across temperatures and EEFs to test if excess entropy scaling relationships for self-diffusion coefficients can extend to ILs in EEF. Fig. 7 shows the total excess entropy computed using eqn (5)versus temperature and across EEFs for replicate 3 as a representative example. For all three replicates we observe a strong excess entropy scaling relationship between structure and dynamics. The excess entropy scaling test for replicate 1 and 2 can be found in Fig. S10 and S11. The excess entropy contribution from each ion was computed by considering translational and configurational excess with eqn (11) and (13). As shown in Fig. 3 and 5, excess entropy is less negative as the magnitude of EEF increases and as temperature increases. According to excess entropy scaling relationships, a less negative excess entropy corresponds to higher dynamics. When the EEF is lower than 1 × 10⁻² V Å⁻¹, it is estimated that excess entropy values at different EEFs are similar to each other. For 1 × 10⁻² V Å⁻¹ and higher EEFs, the estimated excess entropies are much higher for a given temperature. Overall, the estimated excess entropy follows a scaling relationship where excess entropy and reduced combined self-diffusion coefficients across temperature and EEF collapse into a single relationship. The right panel of Fig. 7 shows the scaling relationship, with values for the constants c₁ = 3.9 × 10⁴ and c₂ = 1.4 × 10⁻¹ obtained from the fit to the equation in the right panel with an R² of 0.96 and a RMSE of 0.25. This is consistent with results for excess entropy scaling relationships obtained for many other systems.^8,20 We also tested the scaling relationship for replicate 1 with excess entropy estimates from g(r) processed with KAMEL-LOBE to mitigate any effects of non-converged structural analysis. Fig. S12 and S13 show that the scaling relationship holds with minimal changes, which demonstrates that the computed g(r) results are well converged.⁵⁰ In addition, to the best of our knowledge, this is the first time that excess entropy scaling relationships have been observed for ILs in EEFs. This is an extremely powerful relationship since it shows that the dynamics of an IL can be projected across many thermodynamic states into one relationship. Effectively, the fitted constants and relationship may serve as a dynamical fingerprint of the IL, mapping its entire dynamical landscape across thermodynamic states. We expect this relationship to be transferred to other ILs. Specifically given the correlation between charge density and the response to EEF, we expect that for ILs with ions of low charge density the response of excess entropy to be weaker than for ILs having ions with high charge density.⁷

Fig. 7 Excess entropy scaling test for IL self-diffusion coefficients across temperatures and EEF for replicate 3 as a representative example. (Left) Partial molar configurational and translational excess entropy in [J mol⁻¹ K⁻¹] of the IL versus temperature across EEF. As the temperature and the EEF increase, the magnitude of excess entropy decreases. (Right) Reduced combined self-diffusion coefficient versus partial molar excess entropy. As EEF increases, the magnitude of partial molar excess entropy is less negative which indicates an increase in dynamics. The computed excess entropy and self-diffusion coefficients across EEF and temperatures follow a Dzugutov scaling relationship (black fitted dashed line) with c₁ = 3.9 × 10⁴ and c₂ = 1.4 × 10⁻¹ [J mol⁻¹ K⁻¹]⁻¹ obtained from the fit to the equation in the right panel. The dashed black line is the obtained fit of the equation in the right panel.

Conclusions

Excess entropy scaling relationships can describe the relationship between the dynamics and structure of ILs in EEF and vice versa. Through the estimation of translational and configurational excess entropy, we described and quantified the effect of EEF on the IL structure and showed how structure and dynamics are coupled across all EEFs and temperatures. Through translational excess entropy, we find that the EEF perturbs the IL structure mainly through the [Cl] anion with the [Cl]–[Cl] correlation responding more to the EEF compared to the [C₂C₁im]–[C₂C₁im] correlations, along with effects on the [C₂C₁im]–[Cl] correlation. Translational excess entropy shows that the [Cl]–[Cl] network is affected across short and medium length scales below 10 Å, while the [C₂C₁im]–[C₂C₁im] network is mostly affected beyond 12 Å. In addition, we described and quantified the effect of EEF on IL configurational ordering by estimating configurational excess entropy through the construction of g_conf(r|ω₁) from spatial distribution functions. We observe that at short distances around 4.5 Å, [Cl] localizes in the plane of the imidazolium ring of the [C₂C₁im] cation while other [C₂C₁im] ions localize normal to the imidazolium rings of other [C₂C₁im] cations. We observe that an EEF alters the liquid structure by effectively reducing translational and configurational correlations which are observed with an increase in dynamics. We find that excess scaling relationships describe the increase in dynamics well. The IL structure is mainly described through configurational ordering rather than translational ordering given that configurational excess entropy contributions surpass translational excess entropy. We also observed that the spatial contributions to translational excess entropy are only up to 15 Å, while contributions to the configurational excess entropy go beyond 28 Å. This suggests that IL configurational ordering extends to longer length scales compared to translational ordering. We show that the dynamics of a simple IL with a [Cl] anion can be described with translational and configurational excess entropy of the IL across many temperatures and EEFs. This allows us to obtain material specific constants such as c₁ = 5.0 × 10⁴ and c₂ = 1.5 × 10⁻¹ for [C₂C₁im][Cl] which can serve as a fingerprint that describes the dynamic nature of [C₂C₁im][Cl] across many thermodynamic states. The result of this test is an indication that excess entropy scaling relationships may be able to describe the dynamic properties of other ILs across many thermodynamic states including EEFs.

Author contributions

F. J. C. E. and Y. J. C. conceptualized the research objectives. Data curation was carried out by F. J. C. E. Formal analysis was performed by F. J. C. E., Y. Z., and K. D. Funding acquisition was secured by Y. J. C. and E. M. The investigation was conducted by F. J. C. E. and K. D. Methodology development was led by F. J. C. E. and Y. Z. Project administration responsibilities were managed by Y. J. C., E. M. and F. J. C. E. Resources were provided by Y. J. C., E. M., and Y. Z. Software for determining forcefield parameters was developed and tested by Y. Z. Supervision of the research activities was provided by Y. J. C., E. M., and Y. Z. Validation of results was undertaken by F. J. C. E. and Y. Z. Visualization and preparation of the original manuscript draft were completed by F. J. C. E., while writing—review and editing were collaboratively performed by F. J. C. E., Y. J. C., E. M., Y. Z., and K. D.

Conflicts of interest

The authors have no conflicts to disclose.

Data availability

The data analysis scripts, LAMMPS input files, TRAVIS input files, system data files and forcefield parameters used for this article are available in the Zenodo repository at https://doi.org/10.5281/zenodo.16275259.

Supplementary information (SI) about excess entropy computations for water and results across IL replicates. See DOI: https://doi.org/10.1039/d5cp02782a.

Acknowledgements

This work was supported by the U.S. Department of Energy (DOE) via Subcontract No. 630340 from the Los Alamos National Laboratory, Materials and Chemical Sciences Division, the Department of Education via Grant No. P200A210048, and from the University of Notre Dame (College of Engineering and the Graduate School). The Notre Dame Center for Research Computing provided computational resources. F. J. C. E. would also like to acknowledge Yan Saltar for their contribution in validating the estimation of configurational excess entropy.

References

1. R. Eggert, C. Wadia, C. Anderson, D. Bauer, F. Fields, L. Meinert and P. Taylor, Annu. Rev. Environ. Resour., 2016, 41, 199–222.
2. J. F. Brennecke, J. L. Anderson, G. Belfort, A. Clark, B. Kolthammer, B. Moyer, S. V. Olesik, K. Rosso, M. B. Shiflett, D. Sholl, Z. P. Smith, L. Soderholm, M. Tsapatsis and M. J. Wirth, A Research Agenda for Transforming Separation Science, National Academies of Sciences, Engineering, and Medicine, Washington, DC, 2019.
3. J. F. Brennecke and E. J. Maginn, Ionic liquids: Innovative fluids for chemical processing, 2001, vol. 47.
4. Y. Guan, R. Clark, F. Philippi, X. Zhang and T. Welton, J. Chem. Phys., 2022, 156, 204312.
5. G. Sang and G. Ren, J. Mol. Model., 2018, 24, 1–7.
6. N. J. English, D. A. Mooney and S. O’Brien, Mol. Phys., 2011, 109(4), 625–638.
7. F. J. Carmona Esteva, Y. Zhang, A. S. Arrufat Roman, E. J. Maginn and Y. J. Colón, J. Phys. Chem. B, 2025, 129(28), 7293–7300.
8. J. C. Dyre, Perspective: Excess-entropy scaling, J. Chem. Phys., 2018, 149(21), 210901.
9. X. Ma, W. Chen, Z. Wang, Y. Peng, Y. Han and P. Tong, Phys. Rev. Lett., 2013, 110, 078302.
10. R. Chopra, T. M. Truskett and J. R. Errington, J. Chem. Phys., 2010, 133, 104506.
11. R. V. Vaz, A. L. Magalhães, D. L. A. Fernandes and C. M. Silva, Chem. Eng. Sci., 2012, 79, 153–162.
12. E. H. Abramson, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 062201.
13. Y. Rosenfeld, Phys. Rev. A: At., Mol., Opt. Phys., 1977, 15, 2545–2549.
14. M. Dzugutov, Nature, 1996, 381, 137–139.
15. D. Dhabal, C. Chakravarty, V. Molinero and H. K. Kashyap, J. Chem. Phys., 2016, 145, 214502.
16. A. A. Bertolazzo, A. Kumar, C. Chakravarty and V. Molinero, J. Phys. Chem. B, 2015, 120, 1649–1659.
17. P. Gallo, K. Amann-Winkel, C. A. Angell, M. A. Anisimov, F. Caupin, C. Chakravarty, E. Lascaris, T. Loerting, A. Z. Panagiotopoulos, J. Russo, J. A. Sellberg, H. E. Stanley, H. Tanaka, C. Vega, L. Xu and L. G. M. Pettersson, Chem. Rev., 2016, 116, 7463–7500.
18. A. Barros De Oliveira, E. Salcedo, C. Chakravarty and M. C. Barbosa, J. Chem. Phys., 2010, 132, 234509.
19. O. Lötgering-Lin and J. Gross, Ind. Eng. Chem. Res., 2015, 54, 7942–7952.
20. S. A. Ghaffarizadeh and G. J. Wang, J. Phys. Chem. Lett., 2022, 13, 4949–4954.
21. M. Malvaldi and C. Chiappe, Excess entropy scaling of diffusion in room-temperature ionic liquids, J. Chem. Phys., 2010, 132, 244502.
22. H. Liu and E. Maginn, J. Chem. Phys., 2011, 135, 124507.
23. M. G. Del Pópolo and G. A. Voth, J. Phys. Chem. B, 2004, 108, 1744–1752.
24. M. Agarwal, M. Singh, B. Shadrack Jabes and C. Chakravarty, J. Chem. Phys., 2011, 014502.
25. F. J. Carmona Esteva, Y. Zhang, Y. J. Colón and E. J. Maginn, J. Phys. Chem. B, 2023, 127, 4623–4632.
26. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al. , Comput. Phys. Commun., 2022, 271, 108171.
27. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
28. J. Wang, W. Wang, P. A. Kollman and D. A. Case, J. Mol. Graphics Modell., 2006, 25, 247–260.
29. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, M. B. Frisch, et al., Gaussian 09, Revision A.02, Gaussian, Inc., 2009.
30. C. Schröder, Phys. Chem. Chem. Phys., 2012, 14, 3089–3102.
31. Y. Zhang and E. J. Maginn, J. Phys. Chem. B, 2012, 116, 10036–10048.
32. L. Martinez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
33. M. Součková, J. Klomfar and J. Pátek, Fluid Phase Equilib., 2017, 454, 43–56.
34. K. Jiang, X. Liu, F. Huo, K. Dong, X. Zhang and X. Yao, J. Mol. Liq., 2018, 271, 550–556.
35. M. T. Humbert, Y. Zhang and E. J. Maginn, Mol. Syst. Des. Eng., 2017, 2, 293–300.
36. Y. Zhang, A. Otani and E. J. Maginn, J. Chem. Theory Comput., 2015, 11, 3537–3546.
37. B. Doherty, X. Zhong, S. Gathiaka, B. Li and O. Acevedo, J. Chem. Theory Comput., 2017, 13, 6131–6135.
38. N. Wang and E. J. Maginn, J. Phys. Chem. B, 2024, 128, 871–881.
39. R. Shi and Y. Wang, J. Phys. Chem. B, 2013, 117, 5102–5112.
40. R. Clark, M. Von Domaros, A. J. S. Mcintosh, A. Luzar, B. Kirchner and T. Welton, J. Chem. Phys., 2019, 151, 164503.
41. J. Petravic and J. Delhommelle, J. Chem. Phys., 2003, 118, 7477.
42. N. J. English and C. J. Waldron, Phys. Chem. Chem. Phys., 2015, 17, 12407–12440.
43. B. Lan and Y. Wang, J. Chem. Phys., 2025, 162, 24106.
44. M. T. Humbert, Y. Zhang and E. J. Maginn, J. Chem. Inf. Model., 2019, 59, 1301–1305.
45. M. Brehm, M. Thomas, S. Gehrke and B. Kirchner, J. Chem. Phys., 2020, 152 DOI:10.1063/5.0005078/198859.
46. E. J. Maginn, R. A. Messerly, D. J. Carlson, D. R. Roe and J. R. Elliot, Living J. Comput. Mol. Sci., 2019, 1, 6324.
47. J. Zielkiewicz, J. Phys. Chem. B, 2008, 112, 7810–7815.
48. T. Lazaridis and M. Karplus, J. Chem. Phys., 1996, 105, 4294–4316.
49. J. C. Araque, J. J. Hettige and C. J. Margulis, J. Phys. Chem. B, 2015, 119, 12727–12740.
50. S. A. Ghaffarizadeh and G. J. Wang, J. Chem. Phys., 2023, 224112.

---