The low melting point of room temperature ionic liquids is usually explained in terms of the presence of bulky, low-symmetry, and flexible ions, with the first two factors related to the lattice energy while an entropic effect is attributed to the latter. By means of molecular dynamics simulations, the melting points of 1-ethyl-3-methyl-imidazolium hexafluorophosphate and 1-decyl-3-methyl-imidazolium hexafluorophosphate were determined, and the effect of the molecular flexibility over the melting point was explicitly computed by restraining the rotation of dihedral angles in both the solid and the liquid phases. The rotational flexibility over the bond between the ring and the alkyl chain affects the relative ordering of the anions around the cations and results in substantial effects over both the enthalpy and the entropy of melting. For the other dihedral angles of the alkyl group, the contributions are predominantly entropic and an alternating behavior was found. The flexibility of some dihedral angles has negligible effects on the melting point, while others can lead to differences in the melting point as large as 20 K. This alternating behavior is rationalized by the different probabilities of conformation defects in the crystal.

Ionic compounds are described in textbooks as having high melting points, existing as crystalline solids at room temperature. However, especially after 1990, several ionic compounds with melting points lower than 100 °C or even below room temperature were discovered and were called ionic liquids (ILs). In addition to displaying a low melting point, ILs typically have very low vapor pressure and good conductivity at ambient conditions, characteristics that are desirable for applications like solvents, electrolytes in batteries, lubricants, and heat transfer fluids.1–3 

Low melting points in ILs are usually achieved by the presence of bulky, low-symmetry, and flexible ions. Bulky ions reduce the lattice energy due to the reduction in the electrostatic interaction between oppositely charged ions, thus decreasing the melting enthalpy. Low-symmetry ions also result in less efficient packing in the crystalline phase, contributing to further reduction of the lattice energy. On the other hand, the presence of flexible ions is expected to reduce the melting point due to the increase in the melting entropy, since different conformers may exist in the liquid phase while only some conformations may be allowed in the solid.1,4 However, even with those general guidelines, it is a challenging task to predict the melting point of a new IL or even whether it will be liquid at room temperature.1,4 The presence and the strength of hydrogen bonds between the cations and anions, for instance, affect both the melting point and the viscosity of ILs and were investigated by both quantum and classical calculations.5,6 In the present work, we focused on the question of how the flexibility of alkyl groups in imidazolium-based ILs affects the melting point. A brief review of some known effects of the flexibility over the crystallization, melting, and other properties of ILs and other systems with amphiphilic ions is presented below.

An effect that is attributed to the conformational entropy of alkyl chains in ionic liquids is the so-called “V-shape” seen in a plot of the melting point as a function of the number of carbon atoms, nC, in the linear aliphatic chain of the cations. The melting point initially decreases with the increase in nC until a critical number after which there is a change of tendency and the melting point starts to increase with nC.7–9 The initial decrease is attributed to a dominant effect of the increase in conformational entropy in the liquid phase as the chain gets longer, while the increase for longer chains is explained as a predominance of the enthalpy resulting from van der Waals and hydrophobic interactions between chains.8,9

The possibility of different conformers results in the existence of different polymorphs in the solid phase that differ from each other mainly due to the conformation of some dihedral angles. The IL 1-ethyl-3-methyl-imidazolium dicyanamide, [C2C1IM][DCA], has two solid phases that differ from each other due to the in-plane or out-of-plane conformation of the cation ethyl group, as shown by Raman spectroscopy.10 Computer simulations demonstrate that the conformation of the cation significantly changes the local spatial distribution of the anion.10 The IL 1-butyl-3-methyl-imidazolium chloride, [C4C1IM][Cl], has a monoclinic crystalline form with a gauche dihedral between the first and the second carbon of the butyl group and an orthorhombic form with the same dihedral in the trans conformation.11–13 It was suggested that the dynamics of conformational change in this liquid and in the corresponding bromide, [C4C1IM][Br], leads to supercooling and wide premelting.14,15 Also, different water contents can change the preferred conformation of the cation in the crystal, with the conversion between different polymorphs taking several days in some cases.16 For 1-butyl-3-methyl-imidazolium hexafluorophosphate, [C4C1IM][PF6], three different polymorphs were characterized at different temperatures and pressures. Major differences include the conformation of the cations, the disorder of the anion, and the packing density.17 The protic ionic liquid propylammonium nitrate, [C3H7NH3][NO3], presents the cation in a trans conformation at low temperatures but exhibits a transition to a crystal with gauche conformation upon heating.18 Conformational freedom of ethyl chains was demonstrated to also be an important step during the pre-melting of triethylsulfonium bis(trifluoromethanesulfonyl)imide, [S222][NTf2], upon heating.19 

The flexibility of the anion also affects the melting point if it has multiple conformations. For example, bis(trifluoromethanesulfonyl)imide, [NTf2], is a bulky anion that has two conformers, the trans with point symmetry C2 and the cis with symmetry C1. [NTf2] is known to produce ionic liquids with low melting points.1,4 It has been found that a conformational disordering of the anion takes place before melting, suggesting that the flexibility of the anion contributes to the low melting point.20 Computer simulations showed different free energy surfaces for the two conformers of this anion around the [C2C1IM] cation.21 

The introduction of ether functionalization in cation alkyl chains or perfluoralkyl chains enhances the flexibility of the chain, resulting in both lower viscosity22,23 and lower melting points23,24 than the pure alkyl counterparts.

If, on the one hand, the flexibility affects the physical properties, the chemical environment can also change the conformation or the relative distribution of conformers. The [NTf2] anion, for instance, forms crystals with different conformations depending on the cation.25 Molecular simulations also showed that the relative population of conformers of amphiphilic ions in solution changes due to the differences in the chemical environment. For instance, introduction of urea to an aqueous solution26 or the transference of ions from a hydrophobic to a hydrophilic environment affects the population of different conformers.27 

This brief review shows that the flexibility of both cations and anions is important for the comprehension of the melting point and possible crystalline structures of ILs. However, several questions remain. How much does the flexibility of a given dihedral angle affect the melting point? Or, in other words, if a given flexible dihedral angle is turned rigid, in the sense that it is trapped in a defined conformation, how much will the free energy, the entropy, and the enthalpy difference between the liquid and the solid change? Does this effect depend on the position of the dihedral angle in a given alkyl chain? In the real world, one cannot control the flexibility of a dihedral angle without making chemical modifications like the ether functionalization discussed earlier or the introduction of double bonds. These modifications will introduce other effects alongside with the flexibility change like modifications in the dipole moment or bond lengths. On the other hand, in computer simulation it is possible to restrain a dihedral angle without changing the chemical composition in order to explore a pure flexibility effect. In the present work, we aim to answer these questions by the study of 1-ethyl-3-methyl-imidazolium hexafluorophosphate, [C2C1IM][PF6], and 1-decyl-3-methyl-imidazolium hexafluorophosphate, [C10C1IM][PF6], which will be referred to as [C2] and [C10], respectively, in the rest of this paper. First, the melting point of both ILs was determined using the pseudosupercritical path (PSCP) method.28,29 Then, the effect of flexibility was studied by restraining each dihedral angle one at a time and computing the free energy difference between the solid with and without restraint and between the liquid with and without restraint, which enables a direct calculation of the change of the melting point upon the removal of the flexibility of each dihedral angle.

All the molecular dynamics simulations were performed using the LAMMPS package30 with GAFF (General AMBER force field)31 parameters and the same partial charges employed in previous works.9,32Figure 1 (top) gives the atom names that will be used in the rest of this paper. Note the difference in notation between atom C2 of the cation and system [C2], and the same for C10 and [C10]. The only parameters that were changed from previous works are the parameters for the CR2–N1–C1–C2 dihedral, which were found to give an incorrect description of the secondary minimum as well as an underestimated rotation barrier at 180°. It was reparametrized to fit the results of a potential energy scan done at the MP2/Def2-TVZPD level33 performed with Orca 4.0.1.2.34 The description of the reparameterization and the new parameters are given in the supplementary material. Due to this change, the results presented here differ slightly from the ones of previous works.9,32 The unit cell for both [C2] and [C10] crystals35 is composed of four ion pairs (Fig. 1 bottom). To generate the simulation box for the solid phase, the unit cell of [C2] was replicated five times in the Y and Z directions and four times in the X direction while the unit cell of [C10] was replicated four times in each direction. This resulted in model systems with 400 and 256 ion pairs for [C2] and [C10], respectively. The initial structure for the liquid was prepared by heating the crystal rapidly from 250 K to 600 K at a constant pressure of 1 atm, which ensures spontaneous melting for both compounds. The systems were cooled back to 300 K, and the simulation boxes were deformed to a cube.

FIG. 1.

Top: Structure of the 1-decyl-3-methy-limidazolium cation, [C10], with the atom names labeled. The same names were employed for the 1-ethyl-3-methyl-midazolium cation, [C2]. Bottom: Unit cells used to create the initial structure for the [C2] (left) and [C10] (right) crystals.

FIG. 1.

Top: Structure of the 1-decyl-3-methy-limidazolium cation, [C10], with the atom names labeled. The same names were employed for the 1-ethyl-3-methyl-midazolium cation, [C2]. Bottom: Unit cells used to create the initial structure for the [C2] (left) and [C10] (right) crystals.

Close modal

A cutoff radius of 12 Å was used for both Coulomb and Lennard-Jones interactions as well as for the tethering potential employed in the PSCP method, as described in Sec. II B. Ewald summation was used for long-range electrostatic interactions with a k-space resolution of 0.0001, and a tail correction was employed to account for long-range dispersion interactions. The simulations for the PSCP and the thermodynamic cycle to compute the effect of dihedral angle restraints were performed in the NVT ensemble while the simulations employed for the temperature scans were performed in the NPT ensemble with isotropic coupling for the liquid phases and anisotropic coupling for the solid phases. The Nosé–Hoover thermostat36 with a damping parameter of 100.0 ps was employed in all the simulations, while the barostat with p = 1.0 atm and damping parameter of 100.0 ps was used in the NPT simulations. A time step of 1.0 fs was employed in all simulations.

Structural analyses presented in this paper were performed with the TRAVIS37 software, and the graphical representations were rendered using VMD 1.9.3.38 

The melting point would be largely overestimated if one simply tried to increase progressively the temperature of a bulk crystal in a molecular dynamics simulation. Direct phase coexistence methods are not reliable for computing the melting point of ionic liquids, due to the slow dynamics and complex reordering that occurs in going from the liquid to the crystal. Therefore, special methods are required to compute the melting point. In this work, we employed the PSCP method9,11,28,29,32,39 to compute melting points.

The PSCP method is comprised of four steps to obtain the free energy difference between the solid and liquid phases at a given temperature. These are schematically represented in Fig. 2. In step 1, thermodynamic integration is employed to weaken the non-bonded interactions of the crystal C,

U1=10.9λUvdW+10.9λ2Uelect+λUteth+Ubonded,0λ1.
(1)

At the same time, a tethering potential is turned on to keep the ions vibrating around their lattice points, a state called “weak crystal” (WC),

Uteth=ijaij expbijrij2.
(2)

Notice that the van der Waals interactions, UvdW, are reduced to 10% of their original value while the Coulomb interactions, Uelect, are reduced to 1% of their original value when going from the C (λ = 0) to the WC (λ = 1) state. The potential for bonded interactions, Ubonded, was not changed along any of the thermodynamic integration steps. The tethering potential, Uteth, was applied over the C and N atoms of the cations and the P atoms of the anions. The parameter aij was set to 16.0254 kJ/mol for the cation atoms and 14.0789 kJ/mol for the anion atoms, while bij was set as 0.9 Å−2 for every atom. It was shown in previous works that the choice of the tethering parameters do not affect significantly the computed free energy at the end of the PSCP cycle29 as long as they are strong enough to hold the crystalline structure while reducing the strength of the non-bonded interactions.

FIG. 2.

Illustration of the PSCP method to compute the melting point.

FIG. 2.

Illustration of the PSCP method to compute the melting point.

Close modal

The Helmholtz free energy difference between the C and the WC state was computed by numerical integration of the derivative of the system energy U with the coupling parameter λ using the trapezoidal method as follows:

ΔA=01Uλdλ.
(3)

In step 2, the tethering potential is turned off gradually using the following expression:

U2=0.1UvdW+0.01Uelect+1λUteth+Ubonded,0λ1.
(4)

It is during this step that the melting of the crystal happens. The amorphous state obtained at the end of this step is called “weak dense fluid” (WDF). The free energy change for step 2 is computed again by the integral of U/λ [Eq. (3)].

In step 3, the volume of the system is changed from the crystal density to the liquid density while holding the intermolecular interactions at the same level as at the end of step 2. The resulting state is called the “weak liquid” (WL). The free energy difference between the WDF and WL states is computed via

ΔA3=VCVLpdV,
(5)

where VC and VL are the average volumes of the crystal and liquid obtained in NPT simulations, respectively, and p is the measured pressure of the system.

In step 4, the non-bonded interactions are restored to their full strength,

U4=0.1+0.9λUvdW+0.1+0.9λ2Uelect+Ubonded,0λ1,
(6)

and the liquid, L, is obtained. The Helmholtz free energy variation for this step is obtained again by the integral in Eq. (3) using U4 instead of U1.

The melting point corresponds to the temperature at which the difference in the Gibbs free energy between the solid and liquid phases is zero. To compute the Gibbs free energy at the temperature that the PSCP cycle was performed, the Δ(pV) term between the two phases is added to the sum of Helmholtz free energies computed in the PSCP cycle,11 

ΔGref=ΔA1+ΔA2+ΔA3+ΔA4+pVLVC.
(7)

The Gibbs free energy difference as a function of the temperature was computed from the ΔGref obtained at the reference temperature Tref employed in PSCP and the Gibbs–Helmholtz equation,

ΔGTTΔGrefTref=TrefTΔHTT2dT,
(8)

where ΔH(T) is the enthalpy difference between the liquid and solid, which can be obtained directly from the simulation of both phases at different temperatures.

The enthalpy of both phases and for both the [C2] and [C10] increases linearly with temperature until spontaneous melting is observed for the solid at high temperatures (Fig. S2 in the supplementary material). Discarding the data points above the spontaneous melting, the following linear relation was fit for the enthalpy change:

ΔHT=HLTHCT=ΔH0+ΔCpT,
(9)

where the slope is the difference of heat capacity between the two phases, ΔCp. This linear relation for ΔH(T) was employed to solve the Gibbs–Helmholtz equation, resulting in the closed form equation for the free energy difference between the solid and liquid phases as a function of temperature,

ΔGTT=ΔH01T1TrefΔCp lnTTref+ΔGrefTref.
(10)

Simulations of 4 ns were carried out for each λ in steps 1 and 4 and for each volume in step 3, while 6 ns simulations were performed for each λ in step 2 due to the slower convergence for the λ values in which the transition from the structured to the amorphous state was observed. In each case, the first half of the simulations was discarded in order to guarantee equilibration.

To access the effect of molecular flexibility on melting point, rigidity was added by means of an extra restraining potential applied over the dihedral angles, defined by

Urestr=λK1+cosd,
(11)

where 0 ≤ λ ≤ 1 is a coupling parameter, ϕ is the dihedral angle, and d equals 180° minus the target angle chosen as the most probable angle in the distribution of the dihedral angle in the crystal without restraints. This restrained potential is used in addition to the normal force field parameters, but with a stronger force constant, K = 83.7 kJ/mol. This K is strong enough to enable only small oscillations around the minimum of the potential defined in Eq. (11) at the full strength (λ = 1).

In principle, the same PSCP cycle employed for the unrestrained molecules could be employed for the restrained molecules. However, we used a different thermodynamic cycle to obtain the free energy difference between the new liquid and the new crystal with the restrained dihedral angles, which will be called “Restrained Liquid” (RL) and “Restrained Crystal” (RC) (Fig. 3). In this cycle, in step a, the restrained potential is applied progressively over the crystal C by changing the coupling parameter λ from 0 to 1 in Eq. (11), and the free energy variation of this step, ΔA(a), is computed by numerical integration of Eq. (3). This step is done in the NVT ensemble, and the final state (called RC*) has the same volume as the unrestrained crystal. A volume change is then performed in step b to go from the volume of C to that of RC. The corresponding free energy, ΔA(b), is computed by the integral of −pdV.

FIG. 3.

Thermodynamic cycle used to compute the melting free energy in the system with a restrained dihedral. C and L denote the crystal and the liquid without restraints, RC* and RL* represents the crystal and liquid with the restraints at the same volume of the unrestrained systems, and RC and RL represents the phases with restraints at their respective equilibrium volumes.

FIG. 3.

Thermodynamic cycle used to compute the melting free energy in the system with a restrained dihedral. C and L denote the crystal and the liquid without restraints, RC* and RL* represents the crystal and liquid with the restraints at the same volume of the unrestrained systems, and RC and RL represents the phases with restraints at their respective equilibrium volumes.

Close modal

The same procedure described for the solid is applied to convert L to RL in steps c and d. The free energy difference between C and L is obtained from the PSCP method, which completes the cycle. The convergence of the simulations along this new cycle is faster than the PSCP cycle with smaller statistical error, as will be shown in Sec. III. Only 2 ns simulations were used for each λ in steps a and c or for each volume in steps b and d. Another advantage of introducing this new cycle lies in the interpretation of the results. With this new protocol, one can access directly the effect of the restraints over each phase separately, which would not be possible if the PSCP cycle were employed for the restrained phases.

As in the PSCP method, the Δ(pV) term needs to be added to convert the Helmholtz free energies into the Gibbs free energy as follows:

ΔGref,R=ΔAPSCP+ΔAc+ΔAdΔAaΔAb+pVRLVRC.
(12)

The temperature dependence of enthalpy was determined by performing NPT simulations of the RC and RL as in the case of unrestrained systems. Linear relations were found, and Eq. (10) was used to determine the free energy difference between the RL and RC as a function of temperature.

The value of 83.7 kJ/mol for the constant K in the restraint potential is arbitrary. However, the curves ∂U/∂λ for the solid and the liquid converge to the same values before λ = 0.3 for every dihedral angle studied (Figs. S11–S20 in the supplementary material), which suggests that any K value larger than 25.0 kJ/mol would result in the same free energy difference between RC and RL phases and the same melting point. A value smaller than 25.0 kJ/mol, however, could permit the existence of conformation defects that we want to eliminate by the use of the restraint potential.

The statistical error over the free energy values computed during each step of the thermodynamic cycles was estimated by dividing the equilibrium part of the trajectories into four equal size blocks. Free energy change was calculated using one block of data for each simulation, resulting in four different free energy values for each step. The standard deviation between these blocks was taken as the statistical error for the given step. The statistical error in the free energy difference between C and L and between RC and RL were computed in a similar way using only one block for each step along the cycle [Eqs. (7) and (12)], thus obtaining four different results for the free energy, which enables the calculation of the standard deviation between them. Each of the NPT simulations used to determine the enthalpy variation with temperature was also divided into four blocks. The linear regression defined in Eq. (9) was applied to each set of blocks, resulting in four values of ΔCp and ΔH0 for each system. The standard deviation was then calculated. Finally, using four sets of ΔGref, ΔCp, and ΔH0, four melting points were obtained, and the standard deviation was calculated in a similar manner.

In this section, we will discuss the results for the systems without the dihedral angle restraints, which not only serves as a baseline to discuss the flexibility effect but the free energy obtained from the PSCP method for the flexible molecule [Fig. 2 and Eq. (7)] is also a part of the thermodynamic cycle used for the restrained systems [Fig. 3 and Eq. (12)].

Both [C2] and [C10] crystals do not melt spontaneously when performing NPT simulations at the reference temperatures of 380 K and 425 K, respectively, but significant conformational disorder of the ethyl and decyl groups was observed (see discussion in Sec. III B and III C) as well as rotational disorder of the [PF6] anions. Both effects are expected since similar disorder was observed experimentally even below the melting point for the homologous IL [C4C1IM][PF6]17 and other ILs.20 However, a spontaneous structural change took place for [C10] during the initial relaxation. The N1–C1–C2–C3 and C2–C3–C4–C5 dihedrals were reported to exist as a predominant gauche conformation in the crystal33 (Fig. 1, bottom). However, in the simulation, the system quickly evolved to a different structure in which the trans conformation is preferred. To guarantee that this is not due to improper initial equilibration, we tried to relax the crystal structure at 300 K, which is lower than the experimental melting point. We also tried to apply a restraint to hold these dihedral angles at the initial conformation during the first equilibration simulation and then removed the restraints and performed a second equilibration. In both tests, the crystal evolved to a structure in which every dihedral angle presented the trans as the most favorable conformation. This deviation from experimental crystal structure may be because the force field overstabilizes the trans conformers, an effect observed in a previous work for a different IL with a different force field.9 It is also possible for [C10] to take a crystal structure at higher temperatures that is different than that determined at 173 K.33 Despite this divergence, we proceed with the calculations for this model system and focus on the general features of the dihedral angle flexibility of alkyl groups over the melting of ILs.

Figure 4 shows the simulation results for each step along the PSCP cycle for both [C2] (red curves) and [C10] (black curves). The ∂U/∂λ changes sign during step 1 due to the balance between the non-bonded interactions, which become progressively weaker, and the tethering potential, which becomes stronger as λ increases [Eq. (1)] and eventually becomes dominant over the non-bonded interactions. For [C10], the sign change happens at lower λ, and the final values of ∂U/∂λ are more negative than those for [C2]. As a consequence, the resulting contribution of step 1 per ion pair for the free energy is smaller for [C10] than for [C2] (Table I).

FIG. 4.

Derivative of the internal energy U per ion pair with the coupling parameter λ for steps 1, 2, and 4 along the PSCP cycle and pressure as a function of volume per ion pair for step 3.

FIG. 4.

Derivative of the internal energy U per ion pair with the coupling parameter λ for steps 1, 2, and 4 along the PSCP cycle and pressure as a function of volume per ion pair for step 3.

Close modal
TABLE I.

Results from the PSCP method for [C2] at 380 K and [C10] at 425 K. Values are in kJ/mol per ion pair. The subscript indicates the statistical error in the last digit. For example, 60.304 means 60.30 ± 0.04.

System[C2][C10]
ΔA step 1 (C → WC) 134.053 60.304 
ΔA step 2 (WC → WDF) 88.235 194.868 
ΔA step 3 (WDF → WL) −6.5524 −4.361 
ΔA step 4 (WL → L) −216.994 −251.893 
Δ(pV) (C → L) 0.00201 0.00141 
ΔGref (C → L) −1.266 −1.12 
System[C2][C10]
ΔA step 1 (C → WC) 134.053 60.304 
ΔA step 2 (WC → WDF) 88.235 194.868 
ΔA step 3 (WDF → WL) −6.5524 −4.361 
ΔA step 4 (WL → L) −216.994 −251.893 
Δ(pV) (C → L) 0.00201 0.00141 
ΔGref (C → L) −1.266 −1.12 

A crystalline to amorphous transition occurs during step 2, as can be seen by the drop in ∂U/∂λ between λ = 0.7 and λ = 0.8 in Fig. 4. This step is the one that involves the largest statistical errors and presents slower convergence especially in the region of the transition. A large number of simulations with different λ were carried out at the transition region to ensure an accurate numerical integration. Due to the large number of ions subjected to the tethering potential, the contribution of step 2 is larger for [C10] than for [C2] and is positive in both cases (Table I) since attractive interactions are being removed from the system when turning off the tethering potential.

In step 3, the box dimensions are changed from the crystal phase to the liquid state. Notice that the average pressure is huge at this step since the intermolecular interactions are weakened and the systems tend to expand. The change in the molar volume is larger for [C10], but the pressure is larger for [C2] and the overall contribution for the free energy of this step is more negative for [C2]. In step 4, the Coulombic and van der Waals interactions are restored, leading to the liquid state. This results in a large negative free energy change (Table I) that counteracts most of the changes in steps 1 and 2. Notice that the contribution of step 3 is one or two orders of magnitude smaller than the steps that involve changes in interaction potentials, but is not negligible. The overall free energy changes along the PSCP cycle are of the same order of magnitude of step 3.

Finally, Δ(pV) between the crystal and liquid states was calculated, which makes negligible contribution to the total free energy. The free energy change is negative at the reference temperature for both compounds. This shows that melting should be spontaneous at 380 K for [C2] and at 425 K for [C10]. In the NPT simulations, melting of crystal was observed only at higher temperatures (Fig. S2), showing the need to employ special methods like the PSCP to determine the melting point. Using Eqs. (8) and (10) and results from the PSCP calculations, the melting temperatures were found to be 371 K for [C2] and 409 K for [C10], respectively (Fig. 5).

FIG. 5.

Calculated free energy difference between the liquid and crystal phases. The solid lines stand for the systems without dihedral angle restraints, and the dashed lines stand for the ones with restraints over the CR2–N1–C1–C2 dihedral angle. The vertical lines indicate the computed melting points.

FIG. 5.

Calculated free energy difference between the liquid and crystal phases. The solid lines stand for the systems without dihedral angle restraints, and the dashed lines stand for the ones with restraints over the CR2–N1–C1–C2 dihedral angle. The vertical lines indicate the computed melting points.

Close modal

The predicted melting points for both compounds are higher than the experimental values (Table II). Since the focus in this work is to explore the effects of dihedral angle flexibility, we will use these values as a baseline for further discussion and explore the possible reasons for the overestimation in a separate publication.

TABLE II.

Melting points and changes in enthalpy, entropy, heat capacity, and energy components at the melting temperature without and with the restraint over the CR2–N1–C1–C2 dihedral angle. Values are in K for Tm, J K−1 mol−1 for ΔSm, kJ K−1 mol−1 for ΔCp,m, and kJ/mol for enthalpy and energy components. Experimental data for Tm, ΔHm, and ΔSm are given in parentheses.

System[C2] flexible[C10] flexible[C2] restrained[C10] restrained
Tm 3551 (334,8 336,40 339424092 (307,8 308414302 4052 
ΔHm 19.31 (17.7,8 17.94221.13 (19.4824.91 28.31 
ΔSm 51.63 (53.2,8 52.84251.48 (63.3857.92 69.88 
ΔCp,m 0.0072 −0.0344 0.0422 −0.0655 
Electrostatic 10.244 10.614 12.802 12.122 
van der Waals 7.872 6.975 10.544 12.31 
Bonds 0.1064 0.401 0.0251 0.4871 
Angles 0.382 1.994 0.2572 1.941 
Dihedrals 0.531 1.022 0.961 1.322 
System[C2] flexible[C10] flexible[C2] restrained[C10] restrained
Tm 3551 (334,8 336,40 339424092 (307,8 308414302 4052 
ΔHm 19.31 (17.7,8 17.94221.13 (19.4824.91 28.31 
ΔSm 51.63 (53.2,8 52.84251.48 (63.3857.92 69.88 
ΔCp,m 0.0072 −0.0344 0.0422 −0.0655 
Electrostatic 10.244 10.614 12.802 12.122 
van der Waals 7.872 6.975 10.544 12.31 
Bonds 0.1064 0.401 0.0251 0.4871 
Angles 0.382 1.994 0.2572 1.941 
Dihedrals 0.531 1.022 0.961 1.322 

The effect of dihedral angle flexibility on melting point was studied by the introduction of an extra restraint potential on each dihedral angle at a time [Eq. (11)]. The effect of the CR2–N1–C1–C2 dihedral angle is discussed first due to its different distribution compared to other dihedral angles in the alkyl chain. Without restraints, the maximum distribution of the CR2–N1–C1–C2 dihedral angle in the [C2] crystal is located close to ±90° at high temperature and shifts to larger angles at low temperatures (Fig. 6, solid lines stands for the crystal and dashed lines for the liquid). Another maximum occurs at 0° in the liquid phase for both compounds and the solid phase of [C10], as expected due to the dihedral angle potential energy surface (Fig. S1). For [C2], the differences between the liquid and solid phases are remarkable. In the solid, there is no population at 0°, and the distributions at ± 90° get broader as temperature increases; in the liquid phase, there are three populations with significant superposition, which are independent of temperature. For [C10], on the other hand, the distributions for the solid and liquid phases are very similar. Except for the solid at low temperatures, both phases displayed a small maximum at 0°. The difference between the two compounds is likely related to the different role of van der Waals interactions among the aliphatic tails. The decyl groups in [C10] may maintain a similar cluster structure in both phases while the same effect is not present in [C2]. It is known that the “V-shape” feature in the melting point of imidazolium-based ILs is due to the balance between enthalpy and entropy.7 Relatively, the entropic effect dominates the melting of ILs with small alkyl chains, and the enthalpic effect overcomes the entropic effect for long alkyl group ILs. Different roles of the alkyl group over the structure may also be expected for these two compounds.

FIG. 6.

Distribution of the dihedral angle CR2–N1–C1–C2 in [C2] (top) and [C10] (bottom) at different temperatures without additional dihedral restraints. Solid lines are the distribution in the crystal and dashed lines for the liquid.

FIG. 6.

Distribution of the dihedral angle CR2–N1–C1–C2 in [C2] (top) and [C10] (bottom) at different temperatures without additional dihedral restraints. Solid lines are the distribution in the crystal and dashed lines for the liquid.

Close modal

This dihedral angle has two equally probable populations at ∼± 90°. The dihedral angle restraint [Eq. (11)] was applied to hold half of cations at the positive angle positions and the other half at the negative angle positions in both the liquid and crystal phases. The associated free energy change was calculated along the cycle defined in Fig. 3, and the results are summarized in Table III.

TABLE III.

Calculated Helmholtz free energy change along the thermodynamic cycle of Fig. 3 for the restraint over the CR2–N1–C1–C2 dihedral angle. Free energy values are in kJ/mol, and melting points are in K. The subscripts indicate the computed uncertainty by block average in the last digit.

System[C2][C10]
ΔA step a (C → RC*2.591 6.513 
ΔA step b (RC* → RC) 0.0051 −0.051 
ΔA step c (L → RL*6.651 6.313 
ΔA step d (RL* → RL) −0.0092 −0.0844 
Δ(pV) (RC → RL) 0.00132 0.00021 
ΔGref (RC → RL) 2.919 −1.41 
Tm 4302 4052 
System[C2][C10]
ΔA step a (C → RC*2.591 6.513 
ΔA step b (RC* → RC) 0.0051 −0.051 
ΔA step c (L → RL*6.651 6.313 
ΔA step d (RL* → RL) −0.0092 −0.0844 
Δ(pV) (RC → RL) 0.00132 0.00021 
ΔGref (RC → RL) 2.919 −1.41 
Tm 4302 4052 

Introduction of the dihedral angle restraint produces similar free energy changes for both compounds in the liquid phase (step c in Table III). Compared to the liquid phase, the effect of restraining this dihedral angle is much smaller in the solid phase for [C2] and similar for [C10]. Detailed data are given in Fig. S11 for [C2] and Fig. S12 for [C10] in the supplementary material. The free energy changes due to volume change are negligible in both compounds. As a result, the calculated ΔGref for [C10] with the dihedral angle restraint is very similar to that without the restraints (Tables I and III), and the melting point was found to be 4 K lower, close to the statistical precision limit. For [C2], the restraint causes the ΔGref to increase by more than 4.0 kJ/mol and the melting point to increase by 75 K (Table II and Fig. 5). Note that the statistical errors related to the steps along the cycle shown in Fig. 3 (Table III) are smaller than those along the PSCP cycle (Table I). This demonstrates the advantage of using the cycle in Fig. 3 instead of PSCP for the restrained systems.

The different behavior of [C2] and [C10] upon dihedral angle restraining can be understood by computing the melting entropy and enthalpy. The melting enthalpy can be further decomposed into bonded and non-bonded energy contributions (Table II). Upon imposing a dihedral angle restraint, ΔHm increases by 5.6 kJ/mol and 7.3 kJ/mol for [C2] and [C10], respectively, which means that the restraint stabilizes the solid relative to the liquid from an enthalpic point of view. At the same time, ΔSm increases by 6.3 J K−1 mol−1 and 18.4 J K−1 mol−1, respectively, counteracting the enthalpy change and destabilizing the crystal. For [C2], the enthalpy effect is dominant, and the melting point increases. For [C10], due to the larger entropy change, the effects of enthalpy and entropy cancel out each other, such that the restraint has a small effect on the melting point.

Further analysis reveals that the change in the melting enthalpy arises mostly from the solid phase, and the change in the electrostatic component implies that the organization of the ionic parts is affected by the dihedral angle restraint. This is confirmed by the analyses of the radial pair distribution function between atoms of the cation and anion (Fig. 7). In the liquid phase (dashed lines), there is only a small change of the structure limited to the first coordination shell. For the solid phase, the structural change due to the dihedral angle restraint is more intense and propagates to the second and further coordination shells (solid lines).

FIG. 7.

Radial distribution function, g(r), between the N1 atom of the cation and the P atom of the anion for [C2] at 380 K (top) and for [C10] at 425 K (bottom). The insets are zoom-ins of the second and third coordination shells. The distributions for the systems with restraints over the CR2–N1–C1–C2 dihedral angle are given in red and for the systems without restraints in black.

FIG. 7.

Radial distribution function, g(r), between the N1 atom of the cation and the P atom of the anion for [C2] at 380 K (top) and for [C10] at 425 K (bottom). The insets are zoom-ins of the second and third coordination shells. The distributions for the systems with restraints over the CR2–N1–C1–C2 dihedral angle are given in red and for the systems without restraints in black.

Close modal

Based on the dihedral angle distributions (Fig. 6), one could expect the melting entropy of [C2] to decrease instead of increase upon restraining the dihedral angle. This is because the dihedral angle in the crystal already presents a sharper distribution than in the liquid phase, suggesting a larger entropy change in the liquid phase upon dihedral angle restraining than in the solid phase. If the entropy change ΔΔSm,dih is only related to the dihedral angle distribution, ΔΔSm,dih = ΔSdih,L − ΔSdih,C = (Sdih,RLSdih,L) − (Sdih,RSSdih,S) should be negative. However, ΔΔSm,dih is not the only source of entropy change. Results in Fig. 7 suggest that the intermolecular organization increases more in the solid phase than in the liquid phase. Therefore, the change in the melting entropy due to the dihedral angle restraint can be decomposed into two terms [Eq. (13)], ΔΔSm,dih and ΔΔSm,other, where the first is related directly to the dihedral angle distributions and the latter is due to all the other molecular arrangements in the system. ΔΔSm,other cannot be directly computed without a complete sampling of the phase space. The structural analyses suggest it is positive since the increase in ordering in the crystal due to the restraint is more intense than in the liquid phase. The total ΔΔSm is known to be positive, and ΔΔSm,other must be positive and overcomes the negative term ΔΔSm,dih,

ΔΔSm=ΔSm,restrΔSm,unrestr=ΔΔSm,dihed+ΔΔSm,other.
(13)

For [C10], the changes in the radial distribution functions involving the ionic portions are even more pronounced than for [C2] (Fig. 7), suggesting a larger ΔΔSm,other. The dihedral angle distributions for the two phases are similar in this compound (Fig. 6), indicating a less negative ΔΔSm,dih. Therefore, a larger ΔΔSm is observed for [C10] than [C2] upon the dihedral angle restraint. These qualitative observations agree with the changes in the melting entropy, as shown in Table II.

As a summary for this section, the same dihedral angle displayed different effects on the melting point of the two homologous ILs. The changes in this dihedral angle flexibility affect the intermolecular arrangement, resulting in non-trivial enthalpy and entropy contributions to the melting process. In Sec. III C, the effects of other dihedral angles in the [C10] cation alkyl chain will be discussed.

In this section, the effects of restraints of the other dihedral angles in the decyl group of [C10] are discussed. Each of these dihedral angles will be represented by a number as shown in Fig. 8(a). For example, the dihedral angle N1–C1–C2–C3 (see Fig. 1 for atom names) is numbered as 0, the dihedral angle C1–C2–C3–C4 is numbered as 1, etc.

FIG. 8.

(a) Numbering of the dihedral angles of the aliphatic tail of the [C10] cation, 0 is the N1–C1–C2–C3 dihedral, and the next follows the sequence of carbon atoms in Fig. 1. (b) Fraction of gauche conformers for each dihedral angle in the solid phase at different temperatures. (c) Representations of two layers of the [C10] crystal in the X direction at 300 K (left) and 425 K (right). For better visualization, only the nitrogen and carbon atoms are shown, and the molecules were drawn in different colors depending on the position in the Y direction with red molecules closer and blue farther.

FIG. 8.

(a) Numbering of the dihedral angles of the aliphatic tail of the [C10] cation, 0 is the N1–C1–C2–C3 dihedral, and the next follows the sequence of carbon atoms in Fig. 1. (b) Fraction of gauche conformers for each dihedral angle in the solid phase at different temperatures. (c) Representations of two layers of the [C10] crystal in the X direction at 300 K (left) and 425 K (right). For better visualization, only the nitrogen and carbon atoms are shown, and the molecules were drawn in different colors depending on the position in the Y direction with red molecules closer and blue farther.

Close modal

In the absence of any restraint, typical distributions of alkanes are found for these decyl dihedral angles. They all show a global maximum at 180° corresponding to the trans conformation and symmetrical local maxima at ∼± 70° corresponding to the gauche conformation (Figs. S3–S10 in the supplementary material). Even in the crystal phase, gauche populations have been observed and the gauche fraction tends to increase at higher temperatures. In the liquid phase, no significant change was observed in the distributions with increasing temperature except a small broadening. The gauche fraction in the crystal [Fig. 8(b)] displays an interesting alternating behavior. The even-numbered dihedral angles have smaller populations of gauche defects than the odd-numbered ones. The effect of temperature is stronger over the odd-numbered dihedral angles, which presents a larger increase in the gauche population when increasing the crystal temperature between 300 K and 425 K. Beyond this temperature, the temperature effect becomes small and similar for every dihedral angle.

With increasing temperature, one can notice that for the odd-numbered dihedral angles, the dihedral angle distributions in the crystal converge to those of the liquid phase while the long range ordering remains (Figs. S3–S10). This implies that a significant disordering takes place in the apolar region due to the increase in the less favorable conformation. This behavior is demonstrated for a portion of the [C10] crystal in Fig. 8(c), where the increase in gauche defects and less perfect packing of the decyl groups at higher temperature are observed. On the other hand, the even-numbered dihedral angles have smaller gauche fractions than the liquid phase at all the studied temperatures. The existence of conformational disordering in the crystal even before the melting temperature was already observed experimentally for other ILs.14,15,19

The last dihedral angle (dihedral 7) presents the largest fraction of gauche defects in the crystal. Even in the liquid phase, while the other dihedral angles have nearly the same gauche fraction, this dihedral angle has a larger gauche fraction (Fig. S10). This is because the terminal dihedral angle is less restricted than the others. Due to this peculiarity, some deviations will be noticed for dihedral angle 7 in comparison to the other odd-numbered dihedral angles.

Due to the difference in distribution, it is expected that the flexibility of odd- and even-numbered dihedral angles affect the melting point in different ways. To study this effect, we applied the additional potential defined in Eq. (11) to restrain one dihedral angle each time close to 180° (trans conformation). The computed free energy to restrain each dihedral angle is given in Fig. 9 (top) for the crystal (black lines) and liquid (red lines) phases. The Helmholtz free energy to introduce the restraint is almost constant for every dihedral angle in the liquid phase except for dihedral angle 7 and to a smaller extent for dihedral angle 0. An alternating behavior is observed for the solid phase, which correlates with the gauche fraction shown in Fig. 8(b). The even-numbered dihedral angles that have small gauche fractions show small free energy changes when restrained to the trans conformation because they are already partially restrained due to intermolecular interactions. For a similar reason, the less-restricted odd-numbered dihedral angles experience higher free energy variations for the C → RC conversion.

FIG. 9.

Top: Calculated Helmholtz free energy change with restraints at each dihedral angle for the solid and the liquid of C10. Bottom: Gibbs free energy difference between the restrained liquid and the restrained crystal for each dihedral angle. The dotted line in the bottom panel shows the free energy difference between the liquid and the solid without restraints.

FIG. 9.

Top: Calculated Helmholtz free energy change with restraints at each dihedral angle for the solid and the liquid of C10. Bottom: Gibbs free energy difference between the restrained liquid and the restrained crystal for each dihedral angle. The dotted line in the bottom panel shows the free energy difference between the liquid and the solid without restraints.

Close modal

The free energy difference between the restrained crystal and the restrained liquid was computed according to Eq. (12). The results are shown in Fig. 9 (bottom) where a dashed line was included to represent the free energy difference between the liquid and the crystal in the absence of restraints. As expected, an alternating behavior emerges in an opposite way to that observed in ΔA(C → RC). The even-numbered dihedral angles have large free energy differences because the liquid phase is more affected by the restraint than the solid phase, while the odd-numbered dihedral angles have small free energy changes because the two phases are affected to a similar extent.

It can also be seen in Fig. 9 that while the ΔG(RC → RL) are nearly the same for every even-numbered dihedral angle restrained, ΔG(RC → RL) increases for the odd-numbered dihedral angles when moving further away from the imidazolium ring. If the aliphatic chain is longer, it is likely that ΔG(RC → RL) will converge to the same value at some distance from the imidazolium ring for every dihedral angle except the last one, which is an anomaly. This trend seems consistent with that in the gauche fraction shown in Fig. 8.

The alternating behavior observed in ΔG(RC → RL) (Fig. 9) results in a similar odd–even alternating behavior in the calculated melting point (Fig. 10 and Table IV). Therefore, the flexibility of different dihedral angles of the alkyl groups has different effects on the IL melting point. Restraining the dihedral angle that has a large number of defects in the solid phase does not result in appreciable changes in melting temperature, whereas restraining the ones that tend to have low fraction of defects in the crystal phase can increase the melting point by amounts as high as 20 K.

FIG. 10.

Melting point of [C10] with each dihedral angle restrained. The horizontal dotted line shows the melting point in the absence of restraints.

FIG. 10.

Melting point of [C10] with each dihedral angle restrained. The horizontal dotted line shows the melting point in the absence of restraints.

Close modal
TABLE IV.

Melting points, melting enthalpies, and entropies and differences in energy components between the liquid and crystal at the melting point. Values are in K for melting point, kJ/mol for enthalpy and energy components, and J K−1 mol−1 for the entropy. The subscripts represent the uncertainty in the last digit. Errors smaller than 0.005 are represented as zero.

DihedralNone01234567
Tm 4092 4273 4053 4243 4083 4263 4153 4233 4073 
ΔHm 21.13 20.71 21.314 20.51 21.11 20.71 21.509 20.599 21.608 
ΔSm 51.58 48.51 52.73 48.42 51.72 48.41 51.81 48.71 53.14 
Electrost. 10.614 9.904 10.393 10.135 9.715 10.097 10.485 10.593 10.512 
vdW 6.975 7.938 7.633 7.568 7.865 7.666 7.835 7.125 7.381 
Bonds 0.401 0.380 0.460 0.350 0.470 0.390 0.420 0.370 0.400 
Angles 1.994 1.581 1.951 1.530 2.041 1.571 1.891 1.591 2.092 
Dihedrals 1.022 1.000 0.971 0.940 1.091 0.960 1.051 1.030 1.242 
DihedralNone01234567
Tm 4092 4273 4053 4243 4083 4263 4153 4233 4073 
ΔHm 21.13 20.71 21.314 20.51 21.11 20.71 21.509 20.599 21.608 
ΔSm 51.58 48.51 52.73 48.42 51.72 48.41 51.81 48.71 53.14 
Electrost. 10.614 9.904 10.393 10.135 9.715 10.097 10.485 10.593 10.512 
vdW 6.975 7.938 7.633 7.568 7.865 7.666 7.835 7.125 7.381 
Bonds 0.401 0.380 0.460 0.350 0.470 0.390 0.420 0.370 0.400 
Angles 1.994 1.581 1.951 1.530 2.041 1.571 1.891 1.591 2.092 
Dihedrals 1.022 1.000 0.971 0.940 1.091 0.960 1.051 1.030 1.242 

As opposed to the case of the CR2–N1–C1–C2 dihedral angle, in which both the melting enthalpy and the melting entropy changed appreciably with the restraints, for the other dihedral angles of the decyl group the changes in the melting enthalpy when compared to the unrestrained system are small and the changes in melting entropy are the dominant effect. Because the distribution of each dihedral angle is essentially the same in the liquid (except for dihedral angle 7), the reduction of the entropy in the liquid phase is similar regardless of which dihedral angle restrained. However, in the solid phase the loss of entropy is smaller for the dihedral angles with smaller gauche populations, which implies a smaller ΔSm for the systems with restraints over the even-numbered dihedral angles and therefore causes a larger increase in the melting point.

It is worth mentioning that although it does not have a significant effect on the melting point upon dihedral angle restraining, the melting enthalpy also exhibits an alternating behavior. As shown in Table IV, all the non-bonded and bonded contributions to enthalpy are affected by the dihedral angle restraints. Additional comments regarding these small but systematic variations of the energy components are given in the supplementary material.

The correlation between dihedral angles upon restraining was also studied. It was found that the restraint on one dihedral angle affects the gauche fraction of the other dihedral angles in the crystal phase (black lines in Fig. 11) but not in the liquid phase (red lines in Fig. 11). An alternating behavior in the correlation is also observed. The restraint on one dihedral angle does not affect appreciably the distribution of its immediate neighbor dihedral angles, but instead the gauche population of the next neighbor dihedral angles is suppressed. In other words, the presence of gauche defects for a given dihedral angle n will likely increase the probability of finding gauche defects for dihedral angles n − 2 and n + 2.

FIG. 11.

Fraction of dihedral angles with gauche conformation for each restrained system (solid lines). The gauche fraction for the unrestrained system was included as a reference in each panel (dashed lines).

FIG. 11.

Fraction of dihedral angles with gauche conformation for each restrained system (solid lines). The gauche fraction for the unrestrained system was included as a reference in each panel (dashed lines).

Close modal

The odd–even alternating effect has been reported before and is well known for linear alkanes43 and its α and α,ω-derivates,44,45 where the melting point as a function of carbon number, nC, shows a zig–zag behavior instead of a monotonic shape. The increase in Tm is larger when adding one CH2 to a linear alkane with odd nC than that with even nC. This effect is attributed to the more efficient packing in the crystal phase with even nC. Despite the difficulty to synthesize IL with odd number of carbon atoms in the alkyl tail, examples of similar odd–even effects on the melting point46,47 as well as on density in the liquid phase are known.48 For both neutral molecules and ILs, the odd–even effect tends to disappear as the alkyl chain gets longer, consistent with our observation in melting points as shown in Fig. 10. Such observation makes the tuning of the IL melting point by changing the alkyl chain length and flexibility complicated.

The effects of dihedral angle flexibility on the melting point of alkyl-imidazolium-based ILs were studied using molecular dynamics simulations. Two ILs with different alkyl chain lengths in the cation, 1-ethyl-3-methyl-imidazolium hexafluorophosphate and 1-decyl-3-methyl-imidazolium hexafluorophosphate, were considered. The dihedral angle flexibility around the bond between the imidazolium ring and the alkyl group affects the relative packing and interaction between the charged portion of the cation and the anions, resulting in large variations in both the melting enthalpy and entropy. The resulting effect on the melting point is completely different for the two ILs. For the other dihedral angles of the decyl group, there is no significant effect over the ordering of the ionic part and the effect is mostly entropic, although some systematic enthalpic contributions were also noticed. An “odd–even” effect was found regarding these dihedral angles, which depends on how strong the rotational restraints imposed by the packing in the crystal are. The dihedral angles that are strongly restrained in the solid result in a large conformational entropy gain upon melting; thus, their flexibility contributes significantly to the melting point (up to 20 K). On the other hand, if a dihedral angle has considerable rotational freedom in the crystal, the flexibility affects both liquid and crystal phases in a similar way and no significant effect over the melting point is observed upon imposing a restraint.

If one can increase the dihedral angle rigidity or flexibility in the apolar chain of an IL by means of a chemical modification, for example by introducing an ether group or a double bond, the effect on Tm will depend on the position in the chain. Any practical means to change the flexibility, however, will also change the intermolecular interactions, for instance, by introducing a local dipole by the substituent, which will make the situation complex. Nevertheless, our results can serve as an initial guide to select positions to make chemical modifications to control the melting point of the ionic liquid.

The thermodynamic cycle proposed here to study the effect of flexibility is computationally more efficient than the PSCP method. This new cycle also gives more information as the effect in each phase can be analyzed separately. The same strategy can be used to study other effects over the melting point once the melting free energy in a reference system is determined using the PSCP or other methods.

See the supplementary material for details about the reparametrization of the CR2–N1–C1–C2 dihedral angle, additional comments about the enthalpy effects of decyl group dihedral angle restraints, enthalpy variation with temperature, dihedral angle distributions at several temperatures, raw data used in the thermodynamic integration for introduction of restraints, and distributions of N–C–C and C–C–C angles in the [C10] system.

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

K.B. and M.C.C.R. are indebted to FAPESP (Grant Nos. 2017/12063-8, 2019/04785-9, and 2016/21070-5). M.C.C.R. is also indebted to CNPq (Grant No. 301553/2017-3). Y.Z. and E.J.M. acknowledge support from the Air Force Office of Scientific Research under Contract No. AFOSR FA9550-18-1-0321. We also thank the “Laboratório Nacional de Computação Científica (LNCC/MCTI, Brazil)” for the use of the supercomputer SDumont (https://sdumont.lncc.br).

1.
P. M.
Dean
,
J. M.
Pringle
, and
D. R.
MacFarlane
, “
Structural analyses of low melting organic salts: Perspectives on ionic liquids
,”
Phys. Chem. Chem. Phys.
12
,
9144
9153
(
2010
).
2.
H.
Weingärtner
, “
Understanding ionic liquids at the molecular level: Facts, problems and controversies
,”
Angew. Chem., Int. Ed.
47
,
654
670
(
2008
).
3.
F.
Zhou
,
Y.
Liang
, and
W.
Liu
, “
Ionic liquid lubricants: Designed chemistry for engineering applications
,”
Chem. Soc. Rev.
38
,
2590
2599
(
2009
).
4.
I.
Krossing
,
J. M.
Slattery
,
C.
Daguenet
,
P. J.
Dyson
,
A.
Oleinikova
, and
H.
Weingärtner
, “
Why are ionic liquids liquid? A simple explanation based on lattice and solvation energies
,”
J. Am. Chem. Soc.
128
,
13427
13434
(
2006
).
5.
S.
Zahn
,
G.
Bruns
,
J.
Thar
, and
B.
Kirchner
, “
What keeps ionic liquids in flow?
,”
Phys. Chem. Chem. Phys.
10
,
6921
6924
(
2008
).
6.
P. A.
Hunt
, “
Why does a reduction in hydrogen bonding lead to an increase in viscosity for the 1-butyl-2,3-dimethyl-imidazolium-based ionic liquids?
,”
J. Phys. Chem. B
111
,
4844
4853
(
2007
).
7.
A. I. M. C.
Lobo Ferreira
,
A. S. M. C.
Rodrigues
,
M.
Villas
,
E.
Tojo
,
L. P. N.
Rebelo
, and
L. M. N. B. F.
Santos
, “
On the crystallization and glass-forming ability of ionic liquids: Novel insights into their thermal behavior
,”
ACS Sustainable Chem. Eng.
7
,
2989
2997
(
2019
).
8.
P. B. P.
Serra
,
F. M. S.
Ribeiro
,
M. A. A.
Rocha
,
M.
Fulem
,
K.
Ru̇žička
,
J. A. P.
Coutinho
, and
L. M. N. B. F.
Santos
, “
Solid-liquid equilibrium and heat capacity trend in the alkylimidazolium PF6 series
,”
J. Mol. Liq.
248
,
678
687
(
2017
).
9.
Y.
Zhang
and
E. J.
Maginn
, “
Molecular dynamics study of the effect of alkyl chain length on melting points of [CnMIM][PF6] ionic liquids
,”
Phys. Chem. Chem. Phys.
16
,
13489
13499
(
2014
).
10.
K.
Bernardino
,
T. A.
Lima
, and
M. C. C.
Ribeiro
, “
Low-temperature phase transitions of the ionic liquid 1-ethyl-3-methylimidazolium dicyanamide
,”
J. Phys. Chem. B
123
,
9418
9427
(
2019
).
11.
S.
Jayaraman
and
E. J.
Maginn
, “
Computing the melting point and thermodynamic stability of the orthorhombic and monoclinic crystalline polymorphs of the ionic liquid 1-n-butyl-3-methylimidazolium chloride
,”
J. Chem. Phys.
127
,
214504
(
2007
).
12.
J. D.
Holbrey
,
W. M.
Reichert
,
M.
Nieuwenhuyzen
,
S.
Johnson
,
K. R.
Seddon
, and
R. D.
Rogers
, “
Crystal polymorphism in 1-butyl-3-methylimidazolium halides: Supporting ionic liquid formation by inhibition of crystallization
,”
Chem. Commun.
14
,
1636
1637
(
2003
).
13.
R.
Ozawa
,
S.
Hayashi
,
S.
Saha
,
A.
Kobayashi
, and
H.-o.
Hamaguchi
, “
Rotational isomerism and structure of the 1-butyl-3-methylimidazolium cation in the ionic liquid state
,”
Chem. Lett.
32
,
948
949
(
2003
).
14.
K.
Nishikawa
,
S.
Wang
,
H.
Katayanagi
,
S.
Hayashi
,
H.-o.
Hamaguchi
,
Y.
Koga
, and
K.-i.
Tozaki
, “
Melting and freezing behaviors of prototype ionic liquids, 1-butyl-3-methylimidazolium bromide and its chloride, studied by using a nano-watt differential scanning calorimete
,”
J. Phys. Chem. B
111
,
4894
4900
(
2007
).
15.
H.
Okajima
and
H.-o.
Hamaguchi
, “
Unusually long trans/gauche conformational equilibration time during the melting process of BmimCl, a prototype ionic liquid
,”
Chem. Lett.
40
,
1308
1309
(
2011
).
16.
N.
Kotov
,
A.
Šturcová
,
A.
Zhigunov
,
V.
Raus
, and
J.
Dybal
, “
Structural transitions of 1-butyl-3-methylimidazolium chloride/water mixtures studied by Raman and FTIR spectroscopy and WAXS
,”
Cryst. Growth Des.
16
,
1958
1967
(
2016
).
17.
S.
Saouane
,
S. E.
Norman
,
C.
Hardacre
, and
F. P. A.
Fabbiani
, “
Pinning down the solid-state polymorphism of the ionic liquid [bmim][PF6]
,”
Chem. Sci.
4
,
1270
1280
(
2013
).
18.
L. F. O.
Faria
,
T. C.
Penna
, and
M. C. C.
Ribeiro
, “
Raman spectroscopic study of temperature and pressure effects on the ionic liquid propylammonium nitrate
,”
J. Phys. Chem. B
117
,
10905
10912
(
2013
).
19.
T. A.
Lima
,
V. H.
Paschoal
,
L. F. O.
Faria
, and
M. C. C.
Ribeiro
, “
Unraveling the stepwise melting of an ionic liquid
,”
J. Phys. Chem. B
121
,
4650
4655
(
2017
).
20.
W. A.
Henderson
,
M.
Herstedt
,
V. G.
Young
, Jr.
,
S.
Passerini
,
H. C.
De Long
, and
P. C.
Trulove
, “
Plastic phase transitions in N-ethyl-N-methylpyrrolidinium bis(trifluoromethanesulfonyl)imide
,”
Inorg. Chem.
45
,
1412
1414
(
2006
).
21.
K.
Bernardino
,
K.
Goloviznina
,
M. C.
Gomes
,
A. A. H.
Pádua
, and
M. C. C.
Ribeiro
, “
Ion pair free energy surface as a probe of ionic liquid structure
,”
J. Chem. Phys.
152
,
014103
(
2020
).
22.
M. J.
Monteiro
,
F. F.
Camilo
,
M. C. C.
Ribeiro
, and
R. M.
Torresi
, “
Ether-bond-containing ionic liquids and the relevance of the ether bond position to transport properties
,”
J. Phys. Chem. B
114
,
12488
12494
(
2010
).
23.
J.
Zhang
,
S.
Fang
,
L.
Qu
,
Y.
Jin
,
L.
Yang
, and
S.-i.
Hirano
, “
Synthesis, characterization, and properties of ether-functionalized 1,3-dialkylimidazolium ionic liquids
,”
Ind. Eng. Chem. Res.
53
,
16633
16643
(
2014
).
24.
N.
Terasawa
,
S.
Tsuzuki
,
T.
Umecky
,
Y.
Saito
, and
H.
Matsumoto
, “
Alkoxy chains in ionic liquid anions; effect of introducing ether oxygen into perfluoroalkylborate on physical and thermal properties
,”
Chem. Commun.
46
,
1730
1732
(
2010
).
25.
J. D.
Holbrey
,
W. M.
Reichert
, and
R. D.
Rogers
, “
Crystal structures of imidazolium bis(trifluoromethanesulfonyl)imide “ionic liquid” salts: The first organic salt with a cis-TFSI anion conformation
,”
Dalton Trans.
15
,
2267
2271
(
2004
).
26.
A. F.
de Moura
,
K.
Bernardino
,
O. V.
de Oliveira
, and
L. C. G.
Freitas
, “
Solvation of sodium octanoate micelles in concentrated urea solution studied by means of molecular dynamics simulations
,”
J. Phys. Chem. B
115
,
14582
14590
(
2011
).
27.
K.
Bernardino
and
A. F.
de Moura
, “
Aggregation thermodynamics of sodium octanoate micelles studied by means of molecular dynamics simulations
,”
J. Phys. Chem. B
117
,
7324
7334
(
2013
).
28.
D. M.
Eike
,
J. F.
Brennecke
, and
E. J.
Maginn
, “
Toward a robust and general molecular simulation method for computing solid-liquid coexistence
,”
J. Chem. Phys.
122
,
014115
(
2005
).
29.
Y.
Zhang
and
E. J.
Maginn
, “
A comparison of methods for melting point calculations using molecular dynamics simulations
,”
J. Chem. Phys.
136
,
144116
(
2012
).
30.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
31.
J.
Wang
,
R. M.
Wolf
,
J. W.
Caldwell
,
P. A.
Kollman
, and
D. A.
Case
, “
Development and testing of a general amber force field
,”
J. Comput. Chem.
25
,
1157
1174
(
2004
).
32.
Y.
Zhang
and
E. J.
Maginn
, “
The effect of C2 substitution on melting point and liquid phase dynamics of imidazolium based-ionic liquids: Insights from molecular dynamics simulations
,”
Phys. Chem. Chem. Phys.
14
,
12157
12164
(
2012
).
33.
F.
Weigend
and
R.
Ahlrichs
, “
Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy
,”
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
34.
F.
Neese
, “
The ORCA program system
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
73
78
(
2012
).
35.
W. M.
Reichert
,
J. D.
Holbrey
,
R. P.
Swatloski
,
K. E.
Gutowski
,
A. E.
Visser
,
M.
Nieuwenhuyzen
,
K. R.
Seddon
, and
R. D.
Rogers
, “
Solid-state analysis of low-melting 1,3-dialkylimidazolium hexafluorophosphate salts (ionic liquids) by combined x-ray crystallographic and computational analyses
,”
Cryst. Growth Des.
7
,
1106
1114
(
2007
).
36.
S.
Nosé
, “
A unified formulation of the constant temperature molecular-dynamics methods
,”
J. Chem. Phys.
81
,
511
519
(
1984
).
37.
M.
Brehm
and
B.
Kirchner
, “
TRAVIS—A free analyzer and visualizer for Monte Carlo and molecular dynamics trajectories
,”
J. Chem. Inf. Model.
51
(
8
),
2007
2023
(
2011
).
38.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
, “
VMD—Visual molecular dynamics
,”
J. Mol. Graph.
14
(
1
),
33
38
(
1996
).
39.
G.
Grochola
, “
Constrained fluid λ-integration: Constructing a reversible thermodynamic path between the solid and liquid stat
,”
J. Chem. Phys.
120
,
2122
(
2004
).
40.
V.
Strehmel
,
A.
Laschewsky
,
H.
Krudelt
,
H.
Wetzel
, and
E.
Görnitz
, “
Free radical polymerization of methacrylates in ionic liquids
,” in
Ionic Liquids in Polymer Systems
, edited by
C.
Brazel
and
R.
Rogers
, 1st ed. (
American Chemical Society
,
2005
), pp.
20
21
, available at: https://pubs.acs.org/doi/book/10.1021/bk-2005-0913.
41.
F.
Nemoto
,
M.
Kofu
, and
O.
Yamamuro
, “
Thermal and structural studies of imidazolium-based ionic liquids with and without liquid-crystalline phases: The origin of nano-structure
,”
J. Phys. Chem. B
119
,
5028
5034
(
2015
).
42.
U.
Domanska
and
A.
Marciniak
, “
Solubility of 1-alkyl-3-methylimidazolium hexafluorophosphate in hydrocarbons
,”
J. Chem. Eng. Data
48
,
451
456
(
2003
).
43.
K.
Yang
,
Z.
Cai
,
A.
Jaiswal
,
M.
Tyagi
,
J. S.
Moore
, and
Y.
Zhang
, “
Dynamic odd-even effect in liquid n-alkanes near their melting points
,”
Angew. Chem., Int. Ed.
55
,
14090
14095
(
2016
).
44.
E.
Badea
,
G. D.
Gatta
,
D.
D’Angelo
,
B.
Brunetti
, and
Z.
Rečková
, “
Odd–even effect in melting properties of 12 alkane-α,ω-diamides
,”
J. Chem. Thermodyn.
38
,
1546
1552
(
2006
).
45.
K.
Morishige
and
T.
Kato
, “
Chain-length dependence of melting of n-alcohol monolayers adsorbed on graphite: n-Hexanol, n-heptanol, n-octanol, and n-nonanol
,”
J. Chem. Phys.
111
,
7095
7102
(
1999
).
46.
J. M. S. S.
Esperança
,
M.
Tariq
,
A. B.
Pereiro
,
J. M. M.
Araújo
,
K. R.
Seddon
, and
L. P. N.
Rebelo
, “
Anomalous and not-so-common behavior in common ionic liquids and ionic liquid-containing systems
,”
Front. Chem.
7
,
450
(
2019
).
47.
G.
Adamová
,
R. L.
Gardas
,
L. P. N.
Rebelo
,
A. J.
Robertson
, and
K. R.
Seddon
, “
Alkyltrioctylphosphonium chloride ionic liquids: Synthesis and physicochemical properties
,”
Dalton Trans.
40
,
12750
12764
(
2011
).
48.
G.
Adamová
,
J. N.
Canongia Lopes
,
L. P. N.
Rebelo
,
L. M. N. B.
Santos
,
K. R.
Seddon
, and
K.
Shimizu
, “
The alternation effect in ionic liquid homologous series
,”
Phys. Chem. Chem. Phys.
16
,
4033
4038
(
2014
).

Supplementary Material