Origin of low melting point of ionic liquids: dominant role of entropy

Ionic liquids (ILs) are salts with an extremely low melting point. Substantial efforts have been made to address their low melting point from the enthalpic standpoint (i.e. interionic interactions). However, this question is still open. In this study, we report our findings that entropic (large fusion entropy), rather than enthalpic, contributions are primarily responsible for lowering the melting point in many cases, based on a large thermodynamic dataset. We have established a computational protocol using molecular dynamics simulations to decompose fusion entropy into kinetic (translational, rotational, and intramolecular vibrational) and structural (conformational and configurational) terms and successfully applied this approach for two representatives of ILs and NaCl. It is revealed that large structural contribution, particularly configurational entropy in the liquid state, plays a deterministic role in the large fusion entropy and consequently the low melting point of the ILs.


Supplementary Information Text
Details of IL synthesis. 1,3-dimethylimidazolium iodide ([C 1 mim]I): 1-Methylimidazole (0.268 mol, 22.0 g) and a slight excess amount of iodomethane (0.295 mol, 41.9 g) were dissolved in ethyl acetate. The solution was stirred under an inert atmosphere for 1 hour in an ice bath. The precipitate was washed with ethyl acetate five times, and subsequently recrystallized with acetone. A colorless crystal of [C 1 mim]I was obtained via filtration (yield: 94%). The solution was stirred for 2 hours with light shielding. After the reaction, the precipitant (AgI) was removed by filtration, and the solvent was removed by evaporation. The obtained solid was dissolved in dichloromethane, which produces white precipitant (residual AgI). The precipitant and dichloromethane were removed by filtration and evaporation, respectively. This residual byproduct-removing process was repeated until no precipitant was observed. After the final evaporation, the obtained solid was recrystallized with acetonitrile. A colorless crystal was then obtained via filtration (yield: 64%). 1 H-NMR (DMSO-d 6 ): δ (in ppm) = 9.08 (1H, s, NCHN), 7.67 (2H, t, NCHCH), 3.82 (6H, s, (NCH 3 ) 2 ) 13 C-NMR (DMSO-d 6 ): δ (in ppm) = 137.7 (s, NCHN), 123.9 (s, NCHCH), 36.1 (s, NCH 3 ) Water content: 60 ppm I − content: 320 ppm 1,3-dimethylimidazolium acetate ([C 1 mim]CH 3 CO 2 ): By passing acetic acid (0.235 mol, 14.1 g) aqueous solution through a column filled with the ion exchange resin (Amberlite IRN78, hydroxide form, 55 ml), the ion exchange resin of CH 3 CO 2 − form was obtained. Subsequently, [C 1 mim]I (0.031 mol, 6.95 g) dissolved in distilled water was passed through the column to produce [C 1 mim]CH 3 CO 2 aqueous solution. No detectable iodide salt was confirmed in the solution by the AgNO 3 test. Water was removed from the solution by evaporation and subsequent vacuuming. After washing with ethyl acetate, the obtained solid was recrystallized with acetonitrile to give a colorless crystal (yield: 90%). 1,3-dimethylimidazolium trifluoromethanesulfonate ([C 1 mim]CF 3 SO 3 ): 1-Methylimidazole (0.050 mol, 4.11 g) and a slight excess amount of methyl trifluoromethanesulfonate (0.055 mol, 9.03 g) were dissolved in ethyl acetate. The solution was stirred under an inert atmosphere for 1 day in an ice bath. The solution was evaporated to remove ethyl acetate. The obtained solid was washed with diethyl ether five times, and subsequently recrystallized with acetonitrile. A colorless crystal was obtained as the final product (yield: 95% Brief Theoretical background of two-phase thermodynamic (2PT) approach. The theory of 2PT is briefly described as follows (for details, please see the references 1, 2 ). The total kinetic entropy of a molecule in a liquid state can be divided into translational (S tra ), rotational (S rot ), and intramolecular vibrational (S vib ) contributions. (S1) In 2PT, S tra and S rot where diffusive motions are included are considered as a sum of gaseous and solid components.
The estimation of these entropies is based on density of states function g(v) which is the Fourier transform of the velocity autocorrelation function C(t) of a molecule, where k is the Boltzmann constant, T is temperature, and v is the frequency. C(t) is the sum of the mass-weighted velocity autocorrelation function of atoms, where m is the mass of an atom j and N is the total number of atoms of the systems. Same as the entropy, g(v) is divided into translational (g tra (v)), rotational (g rot (v)), and intramolecular vibrational (g vib (v)) components.
The functions g tra (v) and g rot (v) are determined from autocorrelation functions of center-of-mass velocity and angular velocity of the molecule of interest, respectively. g vib (v) is obtained by the deduction of g tra (v) and g rot (v) from the total density of states function. g tra (v) contains gaseous and solid components.
is expressed by employing a hard-sphere model as, where Δ is the dimensionless diffusivity constant and V is the system volume. The density of states at zero frequency g tra (0) can be determined directly from g tra (v) or via diffusion coefficient D of molecule.
(S11)   tra 12 0 mND g kT  Based on the Carnahan-Starling equation of state, the analytical form of the gaseous translational entropy is expressed as, (S12) where h is the Planck constant, y is the hard-sphere packing fraction, and Z is the compressibility. The estimation of solid translational entropy is based on the harmonic oscillator model.
In the harmonic oscillator model, the canonical partition function of translation Q tra was expressed as,  Figure S3). (S19) Except for Δ 3 A, the Helmholtz energy difference in the PSCP cycle is expressed as, where λ is the alchemical variable ranging from 0 to 1 and U is the potential energy. Starting from the "crystal" state, it is first transformed to the "weak crystal" state. In this step, both LJ (U LJ ) and Coulombic (U Coul ) potentials are weakened, and a tether potential (U tether ) emerges.
The tether potential that binds atoms to lattice points has the Gaussian function form. It was applied for both Na + and Cl − of NaCl. For ILs, the C and N atoms of the cation and the P atom of the anion were used for U tether . The constant a of the cation atoms was 16.0254 kJ mol −1 , and that of the anion atom was 14.0789 kJ mol −1 . 4 The value of 90 nm −2 was used for the constant b of every atom. 4 U bonded is the potential for intramolecular bonds, angles, dihedral angles, and improper angles, which are constant during the PSCP cycle. In step 2, the "weak crystal" is transformed into "weak dense fluid" by removing the tether potential The "weak dense fluid" is then transformed into the "weak liquid". In this step, the cell volume is changed from that of crystal (V cry ) to liquid (V liq ) while the potential is not varied. Therefore, the Helmholtz energy difference (Δ 3 A) is, In step 4, intermolecular potentials are retrieved to the original one, corresponding to the transformation from the "weak liquid" to the normal "liquid" states, as The enthalpy differences between the crystal and liquid states obtained from the NPT MD simulations at various temperatures were fitted with a second-degree polynomial function. Then, equation (S27) becomes, where a, b, and c are the fitting constant. With T m value where ΔG = 0 ( Figure S7) and Δ fus H at the same temperature, Δ fus S is obtained based on equation (1) in the main text. Obtained T m , Δ fus H, and Δ fus S are summarized in Table S7 with reported experimental and MD values. A production run of 2 ns was applied for these simulations 4-7 with 0.1 ps data accumulations.
Conformational analyses were performed on the cations in [C 2 mim]PF 6 or [C 4 mim]PF 6 via 20 ns simulations with 2 ps data accumulations in the NVT ensemble. Once a population of each conformer is estimated from MD trajectories, conformational entropy (S confor ) was calculated, where R is the gas constant and p i is the population of the conformer i. The results are summarized in Tables S8 and S9.
Kinetic entropy (S kin ) estimation. The production runs for 2PT were executed with the velocity Verlet algorithm for 1 ns (a simulation time of 250 ps for 1 block) with 2 fs data acquisition in the NVT ensemble. The convergence of the 2PT method is known to fast (typically ca. 20 ps), 2, 11, 12 and we confirmed that 250 ps is long enough for [C 4 mim]PF 6 ( Figure S10). Translational, rotational, and vibrational density of states functions were obtained using the DoSPT program. 13 The results are displayed in Figures S11-S14, and the obtained numerical values are in Tables S10 and S11 Diffusion coefficient (D) estimation. For diffusion coefficient (D) estimations, production runs of 200 ns with 0.1 ps data acquisition were performed in the NVT ensemble. The diffusion coefficients of the ions were calculated from mean square displacement (MSD) with the Einstein's equation ( Figure S15). The D values are listed in Table S12.          Fig. S14. Absolute kinetic entropies of NaCl (1089 K), [C 2 mim]PF 6 (338 K), and [C 4 mim]PF 6 (275 K). The results in the (A) crystal and (B) liquid states were obtained from the MD simulations (Tables S10 and S11) while that in the (C) gas state was estimated by the DFT calculations (Table S5)