The recalcitrance of lignocellulosic biomass poses a major challenge that hinders the economical utilization of biomass for the production of biofuel, plastics, and chemicals. Ionic liquids have become a promising solvent that addresses many issues in both the pretreatment process and the hydrolysis of the glycosidic bond for the deconstruction of cellulosic materials. However, to make the use of ionic liquids economically viable, either the cost of ionic liquids must be reduced, or a less expensive solvent (e.g., water) may be added to reduce the overall amount of ionic liquid used in addition to reducing the viscosity of the binary liquid mixture. In this work, we employ atomistic molecular dynamics simulations to investigate the impact of water dilution on the overall liquid structure and properties of three imidazolium based ionic liquids. It is found that ionic liquid-water mixtures exhibit characteristics that can be grouped into two distinct regions, which are a function of the ionic liquid concentration. The trends observed in each region are found to correlate with the ordering in the local structure of the ionic liquid that arises from the dynamic interactions between the ion pairs. Simulation results suggest that there is a high level of local ordering in the molecular structure at high concentrations of ionic liquids that is driven by the aggregation of the cationic tails and the anion-water interactions. It is found that as the concentration of ionic liquids in the binary mixture is decreased, there is a point at which the competing self and cross interaction energies between the ionic liquid and water shifts away from a cation-anion dominated regime, which results in a significant change in the mixture properties. This break point, which occurs around 75% w/w ionic liquids, corresponds to the point at which water molecules percolate into the ionic liquid network disrupting the ionic liquids’ nanostructure. It is observed that as the cationic alkyl tail length increases, the changes in the binary mixtures’ properties become more pronounced.

Ionic liquids (ILs) have emerged as a class of highly versatile task specific solvents that are capable of facilitating a myriad of processes that include cellulose dissolution, catalysis, fuel cell applications, and synthesis of novel materials.1–6 In our quest towards a more sustainable future, efforts have been concentrated towards replacing non-renewable fossil-fuel based feedstocks with renewable biomass based alternatives for the production of value added chemicals and transportation fuel.7–9 However, the deconstruction of lignocellulosic materials presents a significant challenge to the economic viability of biomass as a versatile renewable feedstock.10–14 Several ILs have been found that favorably solvate lignocellulosic materials with the imidazolium (Im) based cations paired with anions such as chloride or acetate being the most commonly found ones in biomass pretreatment strategies.11,15–20

In a recent study evaluating the viability of ILs, the key cost-drivers that determine the economic feasibility of cellulosic biofuel production at the industrial scale were identified to be cost of IL production, and employing higher biomass loadings.21 Cellulosic biofuel production processes that involve ILs are classified either into “water-wash,” which involves the removal of ILs with water before enzymatic hydrolysis, or a “one-pot” process, which incorporates both pretreatment and enzymatic hydrolysis in a single step process.22 In previous work, we investigated the conformational energetics of glucose and cellobiose in pure ionic liquids to gain molecular-level insights related to the solvation characteristics and use in effective pretreatment strategies for cellulosic biomass.23–25 However, in the context of the economic feasibility of using ionic liquids, it is imperative that we gain a better understanding of IL-water interactions in order to evaluate the energetic and dynamic behavior of these binary mixtures for the purpose of identifying the optimal mixture ratios.26 While there have been some studies probing the behavior of ionic liquid-water binary mixtures,25,27–29 there still remains a significant scope for an increased understanding of mixture behavior over weight fractions that are relevant experimentally and with respect to the systematic evaluation of the alkyl cation tail length in acetate based ImILs. The primary mode of lignocellulose deconstruction by ImILs is via the disruption of H-bonding interactions.30–32 The questions addressed in this study are crucial to gaining a better understanding of the impact of water dilution and cation tail length on the structural reorganization and its possible impacts on the efficacy of IL components in enabling deconstruction of lignocellulose.

In this study, we employ all-atom molecular dynamics (MD) simulations to evaluate IL-water binary mixtures and the resulting structural and dynamical properties such as molecular ordering, densities, diffusivities, and viscosities. Three imidazolium based acetate ILs with differing side chain lengths (i.e., 1-ethyl-3-methyl imidazolium ([EMIM+]), 1-butyl-3-methylimidazolium ([BMIM+]), and 1-octyl-3-methylimidazolium ([OMIM+])) have been evaluated over four different IL loadings (100%, 75%, 50%, and 25% w/w). These binary mixtures have probed the impact of water dilution on the overall mixture’s characteristics, with special attention on the non-linear transition in system behavior that has been observed at 75% IL loading when transitioning from pure ILs to dilute IL mixtures.27,33–36

The MD simulations for this work were performed using Amber12 and were analyzed using AmberTools12.22 To investigate the characteristics of IL-water binary mixtures and the impact of cationic tail lengths on these characteristics, three ionic liquids were studied: 1-ethyl-3-methyl-imidazolium acetate ([EMIM+] [OAc]), 1-butyl-3-methyl-imidazolium acetate ([BMIM+] [OAc]), and 1-octyl-3-methyl-imidazolium acetate ([OMIM+] [OAc]). Figure 1 shows a schematic of the structures of the cations and anions. The systems studied consisted of ILs in water at various concentrations from 25% to 100%, as shown in Table I. The TIP3P model was used for the water molecules and a charge modified General Amber Force Field (GAFF) was used to describe the ionic liquids.37–39 The Packmol program was used to create the starting coordinates for the pure and binary mixture systems, which were subsequently minimized for 500 steps using the conjugate gradient algorithm.40 The minimization routine was then followed by a temperature-pressure cycling protocol that involved alternating simulations in the isobaric-isothermal ensemble (NPT) at 300 K and 1 atm for 500 ps and the canonical ensemble (NVT) at 400 K for 5 ns, which has been shown to work on viscous and highly entangled polymer systems.23,41–44 This cycling protocol was repeated until the system density converged to within 3% of the reported experimental density for the system at 300 K. All the simulations were performed with the hydrogen atoms restrained using the SHAKE algorithm and with a time step of 2 fs.45 After the cycling simulations, the system was further equilibrated for 10 ns in the NVT ensemble at 300 K, followed by a 100 ns production run in the NVT ensemble at 300 K. After the production run, a 20 ns simulation was performed in the microcanonical ensemble (NVE), from which the self-diffusivity coefficients were calculated.

FIG. 1.

Structures of the imidazolium cations and acetate anion. Nitrogen atoms are depicted in blue, carbons in grey, oxygens in red, and hydrogens in white.

FIG. 1.

Structures of the imidazolium cations and acetate anion. Nitrogen atoms are depicted in blue, carbons in grey, oxygens in red, and hydrogens in white.

Close modal
TABLE I.

Summary of the various ionic liquid (IL) system compositions studied and the cation’s radius of gyration.

ILIL weight percentNumber of IL ion pairsNumber of water moleculesRadius of gyration (Å)
[EMIM+] [OAc25 203 5943 2.63 
50 405 4232 2.63 
75 608 2521 2.63 
100 810 2.64 
[BMIM+] [OAc25 174 5914 3.12 
50 348 4175 3.14 
75 522 2435 3.16 
100 695 3.02 
[OMIM+] [OAc25 136 5876 4.21 
50 271 4098 4.23 
75 407 2320 4.34 
100 542 4.37 
ILIL weight percentNumber of IL ion pairsNumber of water moleculesRadius of gyration (Å)
[EMIM+] [OAc25 203 5943 2.63 
50 405 4232 2.63 
75 608 2521 2.63 
100 810 2.64 
[BMIM+] [OAc25 174 5914 3.12 
50 348 4175 3.14 
75 522 2435 3.16 
100 695 3.02 
[OMIM+] [OAc25 136 5876 4.21 
50 271 4098 4.23 
75 407 2320 4.34 
100 542 4.37 

The cluster analysis of the systems was conducted by first using the Visual Molecular Dynamics (VMD) program to generate a series of images based on system configurations taken every 2 ns from the 100 ns production simulation.46 The series of images for each configuration consisted of 3 Å thick cross-sections that propagated across the entire simulation box. The images were then post-processed with Fiji, a feature-added version of ImageJ, to determine the number, volume, and surface area of domains present for each time frame of the various pure and binary mixture systems.47 In addition, the shapes of the water domains were determined by calculating the sphericity as given by

sphericity=36pVdomain2SAdomain33,
(1)

where Vdomain and SAdomain are the volume and surface area of the domain, respectively.

The local structural order of the IL systems was investigated through the use two order parameters, the Heterogeneity Order Parameter (HOP) and the Orientational Correlation Function (OCF). The HOP is defined as48 

HOP=1NSi=1NSj=1NSexprij22σ2,
(2)

where rij is the distance between sites i and j and σ=L/Ns1/3, L is the simulation box length, and Ns is the number of sites for which the HOP is being computed. The HOP value increases when more sites are closer to each other, resulting in a higher value for more aggregated or clustered geometries. The OCF is used to study the orientational order in the IL cations and is defined based on the vector pointing from the head of the cation to the terminal carbon on the cation tail where the head of the cation is represented by the central carbon (i.e., carbon between the two nitrogen atoms) on the cation’s imidazolium ring. The OCF represents the ensemble averaged correlation, plotted as a function of the distance between the center of mass of the cations, and is defined as49 

OCFr=3u(ri)u(rj)212,
(3)

where u(ri) is the unit vector pointing from the cation head to the cation tail of the cation i, and the correlation is averaged over all cation pairs as a function of distance between cation pairs. The OCF ranges from 0 for completely uncorrelated orientations to 1 for completely correlated cation tails (i.e., parallel or antiparallel orientation).

To study the order in water molecules in binary mixtures, two order parameters were utilized, the Orientational Tetrahedral Order parameter (q) and the Translational Tetrahedral Order Parameter (Sk). The q parameter is one of the most widely used order parameters to study the local structure of water50–56 and is sensitive to the angular order of the water molecules. Its value ranges from 0 for an ideal gas to 1 for a regular tetrahedron. It is computed based on the four nearest oxygen neighbors to the molecule in question and is defined as57,58

q=138j=13k=j+14cosψjk+132,
(4)

where ψjk is the angle between the lines joining the jth and kth oxygen atoms among the 4 nearest neighboring atoms to the molecule being considered. The Sk order parameter, although less frequently used in the literature,56,58–60 has been shown to be more sensitive to density fluctuations than q61 and increases when the local tetrahedral order increases, up to a maximum value of 1 for a tetrahedron. The Sk order parameter measures the variance of radial distances for the four nearest neighboring water oxygens to the water molecule being considered and is defined as

Sk=113k=14rkr¯24r¯2,
(5)

where rk is the distance from the water oxygen being considered to the kth neighboring atom, and r¯ is the mean of the distances to all four neighboring water oxygens.

The results of the temperature-pressure cycling protocol yielded densities that were within 3% of experimentally determined densities. In Figure 2, it is evident that the densities calculated from our simulations are in good agreement with available experimental values as well as with previously published simulations for [EMIM+] [OAc] and [BMIM+] [OAc]. It is noted that as of yet, experimental densities are not available for [OMIM+] [OAc], although the trend of decreasing density with increasing cation tail length agrees with previous experimental and theoretical work.27,36 The density as a function of IL content exhibits an unusual parabolic trend (i.e., density turnover) with the densities peaking at an IL concentration of about 75% IL in water. Similar trends of an increase in the density of ILs with addition of water (75%–100% in Figure 2) have been reported in experimental studies33–36 of [RnMIM+] [OAc] IL-water systems, where Rn represents the alkyl tail (n = 2, 3, 4, 5, and 6). Previous work has attributed this parabolic behavior to a change in the molecular-level ordering and ion-water interactions exhibited by the IL ion pairs.27 In order to obtain a more detailed understanding of this unusual behavior of ILs, the effect of water dilution on the local structure and behavior of ILs is quantified using different analyses, which are discussed in greater detail in Secs. III B and III C.

FIG. 2.

Pure IL and water-IL binary mixture densities. The current simulation results for [EMIM+] [OAc] (black square), [BMIM+] [OAc] (red square), and [OMIM+] [OAc] (blue square) are represented with a solid line connecting the data points. Experimental work for [EMIM+] [OAc]32 (black circle) and [BMIM+] [OAc]1–4,31 (red circle) is represented with a dotted line connecting the data points, while previous simulation work25 for [EMIM+] [OAc] (black triangle) and [BMIM+] [OAc] (red triangle) is represented with a dashed line connecting the data points. Error for the depicted data is denoted by the size of the symbols.

FIG. 2.

Pure IL and water-IL binary mixture densities. The current simulation results for [EMIM+] [OAc] (black square), [BMIM+] [OAc] (red square), and [OMIM+] [OAc] (blue square) are represented with a solid line connecting the data points. Experimental work for [EMIM+] [OAc]32 (black circle) and [BMIM+] [OAc]1–4,31 (red circle) is represented with a dotted line connecting the data points, while previous simulation work25 for [EMIM+] [OAc] (black triangle) and [BMIM+] [OAc] (red triangle) is represented with a dashed line connecting the data points. Error for the depicted data is denoted by the size of the symbols.

Close modal

1. Radial distribution functions (RDFs)

To obtain a more detailed understanding of the local ordering structure of the IL-water systems and the interactions between the cations, anions, and water molecules, various radial distribution functions (RDFs) were computed for the components of the system. The coordination numbers for the RDFs were also calculated by integrating the first peak of the RDFs and are shown for all the systems studied. In order to maintain consistency in our comparisons, the coordination numbers were calculated by considering a standard distance across all dilutions. The specific distance cutoffs considered for coordination number calculations are listed in Table I of the supplementary material. In our previous work, we have discussed the RDFs for pure ILs23,24 and therefore will be focusing the current discussion on the binary mixtures and how they differ from the pure IL systems. Figure 3 shows the RDFs for the ion pairs of [OMIM+] [OAc], for the pure IL and the IL-water binary mixtures. Due to the similarity in the general structure of the RDFs for the three ILs investigated, only one set of RDFs is presented. The RDFs for [EMIM+] [OAc] and [BMIM+] [OAc] can be found in the Figure 1 of the supplementary material. From Figure 3, we can see that at 100% IL there are two prominent peaks in the RDFs, which is in good agreement with previous reports in the literature.23,27 This relatively strong bimodal distribution signifies a large amount of local ordering among the cations and anions in the first few solvation shells. Upon addition of water, it is observed that the solvation structure between the cation and anion is weakened as indicated by the reduced amplitude of the first two peaks in the 75%, 50%, and 25% cation-anion RDFs. Hence, the addition of water disrupts the strong cation-anion association and ordering. The differences between the [EMIM+], [BMIM+], and [OMIM+] RDFs suggest that the length of the cation tail has a major impact on the solvation structure as seen in Figure 4. For the [EMIM+] [OAc] and [BMIM+] [OAc], there is a more drastic effect of dilution with respect to the cation-anion peaks, wherein as water is added, the amplitude of the first two peaks is significantly reduced and ultimately results in a single broad distribution. In the [OMIM+] [OAc] IL, the reduction in the peak amplitudes is not as severe and even at the 25% IL mixture two distinct peaks are present. This suggests that the longer cation tails in [OMIM+] [OAc] contribute to the preservation of the local ordering, even at low IL concentrations, due to the larger number of hydrophobic interactions in the cation tail. The 100% IL cation-cation and anion-anion RDFs (Figure 3), on the other hand, are characterized by broad first peaks representing the first solvation shell and a shoulder feature indicating subsequent solvation shells. Upon dilution, the cation-cation RDFs do not change significantly, although the first peak is observed to move to a larger distance, which is expected due to the increasing presence of intervening waters. The cation-cation RDFs converge to a low structure distribution between 50% and 75% IL. The anion-anion RDFs, however, exhibit a very interesting trend. As the IL concentration decreases to 75%, there is a sharp, high intensity peak at 6 Å for all the ILs studied, as well as clear structures indicating subsequent solvation shells. Hence, the anion-anion RDF change from a broad distribution at 100% IL to distinct peaks at 75% indicates that the hydrophilic anions interact strongly with water, hence inducing a local ordering in the anion solvation environment when in the presence of water. This is also evident in the increase in the coordination number for anion-anion interactions from 1.64 to 2.04 in going from 100% to 75% IL concentration. As the concentration of water is further increased, intensity of the first peak as well as the coordination numbers for the anion-anion interactions decreases, hence signifying that the local ordering in the anions is at a maximum at 75% IL concentration. This could be attributed to the fact that at lower water concentrations (75% IL), there is a strong interaction between the water and anions in the system. As the amount of water increases, water-water hydrogen bonding becomes more prevalent, thus leading to a reduction in the anion-water interaction and hence the anion structuring. This is exacerbated by the fact that at low IL concentrations there are significantly fewer numbers of anion-anion interactions as they are increasingly diluted in the aqueous solvent. Hence, the anions in the IL-water mixture show a favorable attraction towards the water, and the anion-anion RDFs converge to an ordered solvation structure in going from 100% to 75% IL concentration. This trend is consistent with our previous results that identified the presence of a strong IL-IL interaction at high concentrations of IL and a subsequent disruption of this interaction at IL concentrations below 75%. Also, combining this observation with the cation-cation RDF observations, it can be concluded that the strong polar interactions between the anion and water at low water concentrations (75% IL) induce not only anion-anion structuring but also impact the cation headgroups, hence causing a structured polar domain at low water concentrations.

FIG. 3.

RDFs for cation-anion (black solid line), cation-cation (red dashed line), and anion-anion (blue dotted line) interactions for (a) 100%, (b) 75%, (c) 50%, and (d) 25% IL dilution for [OMIM+] [OAc]. The coordination numbers are reported in the inset.

FIG. 3.

RDFs for cation-anion (black solid line), cation-cation (red dashed line), and anion-anion (blue dotted line) interactions for (a) 100%, (b) 75%, (c) 50%, and (d) 25% IL dilution for [OMIM+] [OAc]. The coordination numbers are reported in the inset.

Close modal
FIG. 4.

RDFs for cation-anion for [EMIM+] [OAc] (red solid line), [BMIM+] [OAc] (black dashed line), and [OMIM+] [OAc] (blue dotted line) interactions for (a) 100%, (b) 75%, (c) 50%, and (d) 25% IL dilution. The coordination numbers for the first solvation shell are reported in the inset.

FIG. 4.

RDFs for cation-anion for [EMIM+] [OAc] (red solid line), [BMIM+] [OAc] (black dashed line), and [OMIM+] [OAc] (blue dotted line) interactions for (a) 100%, (b) 75%, (c) 50%, and (d) 25% IL dilution. The coordination numbers for the first solvation shell are reported in the inset.

Close modal

To study these ion-water interactions in further detail, the ion-water RDFs were computed, as shown in Figure 5. Again, only the RDFs for [OMIM+] [OAc] are shown, as the RDFs for all the ILs studied showed a similar overall structure (specific RDFs can be found in Figure 2 of the supplementary material). From the cation head-water RDFs, anion-water RDFs, and the corresponding coordination numbers, it is clearly seen that the ion-water interaction is the strongest at 75% and decreases with decreasing IL concentration, as exhibited by the reduction in the RDF peak amplitude when going from 75% to 25% IL concentration. Furthermore, the anion-water peak for all concentrations is higher than the water-water peak, and it increases as the IL concentration is increased, indicating the strong anion-water interactions, which peak at 75% IL. It is further observed that while the water-water coordination number is higher than the anion-water coordination number at 25%, this difference decreases with an increase in IL concentration, and at 75%, the coordination number for anion-water interactions is higher than the water-water coordination number. This further supports the hypothesis that the anion-water interactions increase with increasing IL concentration up to a maximal value at around a concentration of 75% IL. Combined with our previous results for the RDFs, it is concluded that a maximal local structural ordering in the IL-water binary mixtures is reached at around 75% IL. The cation tail-water RDFs show that the cation tail-water interactions are relatively weak, with the intensity of these interactions decreasing with increasing IL concentration, implying cation tail aggregation in the presence of water. The tail aggregation is seen to increase with increasing IL concentration. Furthermore, there is a significant effect of tail length observed for [OMIM+] [OAc] as compared to [EMIM+] [OAc] and [BMIM+] [OAc], where the amount of tail aggregation increases with an increase in tail length, showing that while there may not be a significant amount of tail aggregation for [EMIM+] [OAc], the tails begin to aggregate for [BMIM+] [OAc] and aggregation is maximum for [OMIM+] [OAc] (supplementary material, Figure 3).

FIG. 5.

RDFs for [OMIM+] [OAc] for cation head-water (red dashed line), cation tail-water (green dashed-dotted line), anion-water (black dotted line), and water-water (blue solid line) interactions for (a) 75%, (b) 50%, and (c) 25% IL dilution. The coordination numbers are reported in the inset.

FIG. 5.

RDFs for [OMIM+] [OAc] for cation head-water (red dashed line), cation tail-water (green dashed-dotted line), anion-water (black dotted line), and water-water (blue solid line) interactions for (a) 75%, (b) 50%, and (c) 25% IL dilution. The coordination numbers are reported in the inset.

Close modal

Finally, to understand the interaction between the ions and water in the IL-water solutions, we computed the interaction energies between the ions and water by using the reversible work theorem,

gr=expArkBT,
(6)

where g(r) is the radial distribution function, kB is the Boltzmann constant, T is the temperature, and A(r) is the Helmholtz free energy (i.e., simulations were in the NVT ensemble). The results for these calculations are shown in Figure 6. The interaction energies were calculated from the corresponding RDFs, and a bootstrap error analysis on the RDFs revealed that the RDFs had already converged well, and hence the error on the interaction energies was low (of the order of 0.1 kcal/mol). This error was accepted to be satisfactorily low considering the high number of interactions and the significant amount of time the values were averaged over. Also, for the interaction energy calculation, we used the difference between the maximum of the first peak of the RDFs to the value of the RDF at infinite separation (∼1), which represents the thermodynamic energy required to separate the interacting atoms to an infinite distance. From the figure, we can see that the ion-water interaction energies increase with increasing IL concentration. In addition, the cation-water interaction energies are very low as compared to the anion-water interactions, suggesting that the anions have a much stronger interaction with the water as would be expected. The anion-water interaction energy at low IL concentrations is slightly higher than the water-water interaction energies, with this difference increasing with increasing IL concentration. At 75% IL, we see that the anion-water interactions are significantly stronger than the water-water interactions. This result agrees with our previous RDF results that suggested that anion-water interactions are dominant at 75% IL and decrease with decreasing IL concentration, hence leading to the highest local structural ordering among ions at around 75% IL. This results also explains the density turnover mentioned before, where the density for the ILs shows a maximum at around 75% IL. This can now be attributed to the relatively high anion-water interactions at around 75%, hence leading to high local structural ordering and hence a higher density than that observed at other IL concentrations. Previous work suggested that as the length of the cation tail increased to 6 and 8 carbons, the density turnover moved to a lower IL concentration27,62 (i.e., more water was required to disrupt the IL structure). This phenomenon cannot be observed in our results because we did not have simulations for IL concentrations between 50% and 70%. However, this can be explained from Figure 6 by comparing the anion-water and cation tail-cation tail interaction energies. From the figure, we can see that for [EMIM+], the tail-tail interaction energies are much lower than the anion-water interaction energies. This implies that for smaller tails, the anion-water electrostatic interaction induced structuring of the polar domain is the dominant mechanism for the density turnover. However, as the chain length increases, the tail-tail interaction energies increase. For [OMIM+], the tail-tail interaction energies are comparable to the anion-water energies at 50% and 75% IL concentration, while it is higher than the anion-water energy at 25%. This implies that as the tail length increases, the effect of tail aggregation on the local ordering becomes more significant and ultimately the dominant mechanism. Hence, in IL-water systems, there are two factors that contribute to local ordering in the IL structure, anion-water interactions leading to structuring of the polar domain and cation tail aggregation leading to structuring of the non-polar domain. For longer tails, the effect of the latter is larger, hence causing a competition between a breakdown of the structure imposed by water penetrating the polar domain at low IL concentrations and the strengthening of the non-polar domain due to tail aggregation. This tail aggregation driven non-polar structure formation resists a breakdown in the IL structure up to higher water concentrations for longer tails, hence explaining why the density turnover is observed at lower IL concentration for longer tails.

FIG. 6.

Interaction energies for cation head-water (black square), anion-water (red circle), cation tail-cation tail (green diamond), and water-water (blue triangle) interactions for (a) [EMIM+] [OAc], (b) [BMIM+] [OAc], and (c) [OMIM+] [OAc].

FIG. 6.

Interaction energies for cation head-water (black square), anion-water (red circle), cation tail-cation tail (green diamond), and water-water (blue triangle) interactions for (a) [EMIM+] [OAc], (b) [BMIM+] [OAc], and (c) [OMIM+] [OAc].

Close modal

2. Clustering analysis

In order to obtain further insight into the local nanoscale ordering and how the addition of water to ILs impacts cluster distributions, an analysis was performed on the water clusters for all systems. Figure 7 indicates that there is no significant change in water cluster size from pure water to 50% IL binary mixture for the three ILs investigated (i.e., there exists one continuously connected water cluster). However, at 75% IL it can be seen that there is a noticeable deviation from the more dilute systems. When the system contains at least 75% IL the presence of smaller water domains becomes more prevalent. There is also a trend observed related to the increase in the cation tail size, with an increase in tail length resulting in a smaller deviation from the dilute systems. This is also consistent with our previous observations of an increased tail clustering with increasing tail lengths, which would result in larger cation domains, leading to larger water domains.

FIG. 7.

Size of largest water domains for [EMIM+] [OAc] (black square), [BMIM+] [OAc] (blue circle), and [OMIM+] [OAc] (red triangle) systems. All error bars are smaller than the symbol size.

FIG. 7.

Size of largest water domains for [EMIM+] [OAc] (black square), [BMIM+] [OAc] (blue circle), and [OMIM+] [OAc] (red triangle) systems. All error bars are smaller than the symbol size.

Close modal

In addition to the largest domain size, the shape of the water clusters was also evaluated. Figure 8 shows the normalized sphericity distribution for each of the IL systems studied. It is noted that a sphericity of 1 represents a perfect sphere, while more elongated asymmetrical shapes will have a low sphericity. It can be seen that all systems have a peak at about 0.1, which represents the large water domain in each system as seen in Figure 7. When other peaks appear they are in the range of 0.4–0.5 for all systems, suggesting that the shape of smaller domains for different ionic liquids is similar, and it is only the frequency of their occurrence that tends to change, as indicated by the amplitude of the curves. At the 25% and 50% ionic liquid levels, the curves are all very similar, with few smaller domains, which is consistent with the results depicted in Figure 7. However, at 75% IL, the average number of clusters having a sphericity value of 0.4–0.5 is much higher than that at lower dilutions, and this reduces with increasing tail length. As with the previous conclusion from this section, this trend is seen because the hydrophobic tails of the cations tend to cluster together in these ionic liquid-water mixtures, with longer tails able to do so more effectively.

FIG. 8.

Normalized sphericity distribution of water domains in water-IL systems for 75% (blue dotted lines), 50% (red dashed lines), and 25% (black solid lines) IL for (a) [EMIM+] [OAc], (b) [BMIM+] [OAc], and (c) [OMIM+] [OAc] systems.

FIG. 8.

Normalized sphericity distribution of water domains in water-IL systems for 75% (blue dotted lines), 50% (red dashed lines), and 25% (black solid lines) IL for (a) [EMIM+] [OAc], (b) [BMIM+] [OAc], and (c) [OMIM+] [OAc] systems.

Close modal

3. Order parameter analysis

To further investigate the local structure and the molecular-level interactions between various components in the IL-water systems, we selected two sets of order parameters to study the structure of both the ions and water molecules in the binary mixtures. The first set of order parameters, the Heterogeneity Order Parameter (HOP) and the Orientational Correlation Function (OCF), was used to study the local structure of the ILs. Figure 9 shows the results for the HOP analysis and indicates that for the 100% IL concentration the cation heads and anions have relatively low values indicative of a low ordered homogeneous solvent. The HOP value does not change much for smaller cation tail lengths, but increases slightly for [OMIM+] [OAc] due to the increased cation-cation interactions and increased anion-water interactions. The cation tails, however, show a clear trend of an increasing HOP with an increase in cation tail length. This can be attributed to the increased aggregation of cationic tails, hence leading to a closer proximity of cation tail atoms, resulting in a higher HOP value. This observation of an increased spatial heterogeneity with an increase in tail length in ILs is consistent with previous studies in the literature,49,63,64 which have also attributed the increase in the HOP to an increased tail aggregation with increasing cation tail length. As the concentration of ILs is changed, the cation tail HOP is observed to undergo changes that increase in magnitude as the size of the tail increases. It is observed that [EMIM+] [OAc] does not exhibit a significant trend with an increase in IL concentration, while [BMIM+] [OAc] indicates that the HOP increases with increasing IL concentration up to a concentration of 75% and then decreases with an additional increase in IL concentration. For [OMIM+] [OAc], a similar trend is seen; however, the increase in HOP with concentration is significantly more drastic than for the two other ILs. It is observed that the HOP increases until a value of 75% IL, beyond which the HOP decreases drastically between 75% and 100% IL. This result further corroborates the conclusion that there is a higher ordering in local structure for the IL in the presence of water, with a maximum in ordering being reached at about 75%, arising from increased aggregation of the cation tails.

FIG. 9.

HOP for (a) 100% IL as a function of tail length for anions (black square), cation heads (red circle), and cation tails (blue triangle) and (b) cation tails as a function of concentration for [EMIM+] [OAc] (black square), [BMIM+] [OAc] (red circle), and [OMIM+] [OAc] (blue triangle).

FIG. 9.

HOP for (a) 100% IL as a function of tail length for anions (black square), cation heads (red circle), and cation tails (blue triangle) and (b) cation tails as a function of concentration for [EMIM+] [OAc] (black square), [BMIM+] [OAc] (red circle), and [OMIM+] [OAc] (blue triangle).

Close modal

Figure 10 shows the results for the OCF values computed for the various IL systems studied. At small distances, the values of OCF are seen to be negative. This would imply that the tails are perpendicular to each other and is expected because at small distances, the cation heads will repel each other. In addition, at low distances the tails cannot be anti-parallel to each other as they would prefer to assume a relative orientation that is closer to a perpendicular orientation to reduce cation-cation charged repulsion. The first peak in the OCF values corresponds to the correlation between two closest cations that adopt a parallel or antiparallel orientation, while the smaller peaks at larger distances correspond to a long range order. The figure shows that this long range order is present to a larger extent for ILs with longer cation tails as seen by higher amplitude features at larger distances. Two distinct orientational preferences are observed for [EMIM+] [OAc], with larger OCF values for 25% and 50% IL mixtures and a significant reduction in the OCF for the 50% and 75% IL mixtures. This would indicate that in dilute IL mixtures there is a tendency for the small cation tails to align in a parallel or antiparallel fashion, but this is quickly eliminated upon increasing IL concentrations. It is also noteworthy that for the 25% IL system there appears to be no perpendicular orientation of the cations due to the low concentration and relative freedom to reorient to minimize cation-cation repulsion. In [BMIM+] [OAc] all dilutions studied indicate a reduced ability to align the tails in a parallel or antiparallel orientation as compared to the 100% IL. The trend of higher OCF values at higher IL concentrations is further amplified when evaluating [OMIM+] [OAc]. For the longer cation tail IL, it is observed that the 100% and 75 IL systems possess the largest OCF values with the maximum orientation occurring at 75% IL. The OCF is then observed to decrease upon increasing water content. This result can be explained by the increase in tail aggregation in the presence of water upon increase in IL concentration, again, with a maximum occurring around 75% IL.

FIG. 10.

Orientational Correlation Function (OCF) values for 25% (black solid line), 50% (red dashed line), 75% (blue dotted line), and 100% (green dashed-dotted line) IL concentrations computed for (a) [EMIM+] [OAc], (b) [BMIM+] [OAc], and (c) [OMIM+] [OAc].

FIG. 10.

Orientational Correlation Function (OCF) values for 25% (black solid line), 50% (red dashed line), 75% (blue dotted line), and 100% (green dashed-dotted line) IL concentrations computed for (a) [EMIM+] [OAc], (b) [BMIM+] [OAc], and (c) [OMIM+] [OAc].

Close modal

The second set of order parameters consisted of the Orientational Tetrahedral Order (q) and the Translational Tetrahedral Order (Sk), which were used to investigate the local structure in the water molecules. Figure 11 reports the averages for both these order parameter values for varying IL concentration. The figure shows that on addition of IL to water, the local order of the water is reduced as compared to bulk water as the water favorably interacts with the anion. The value of the order parameters is observed to increase with an increase in cation tail length, which is hypothesized to be due to an increase in cation tail ordering, which disfavors water and anions in the cation-tail rich (i.e., hydrophobic) regions. Figure 12 shows the probability distribution of q for the [OMIM+] [OAc]-water systems compared to the values for bulk water. The distributions for [EMIM+] [OAc] and [BMIM+] [OAc] are similar in structure to that for [OMIM+] [OAc] and are presented in Figure 4 of the supplementary material. From the figure, it is clear that the orientational order of the water molecules decreases on addition of IL, and as the IL concentration is increased the distribution goes from a more water-like distribution at 25% IL to a broader distribution at 75% water. The broadening of the distribution signifies a depletion in the ordered water structures at 75% IL concentration as compared to the lower IL concentrations, again implying a higher IL ordering at 75% IL concentration. As mentioned before, this could be due to the higher anion-water interactions as well as due to the higher cation tail aggregation around 75% IL, hence leading to a reduction in the structural ordering in water.

FIG. 11.

Average values of water order parameter values for IL-water systems studied for [EMIM+] [OAc] (black square), [BMIM+] [OAc] (red circle), and [OMIM+] [OAc] (blue triangle) for (a) Orientational tetrahedral order (q) and (b) Translational tetrahedral order (Sk).

FIG. 11.

Average values of water order parameter values for IL-water systems studied for [EMIM+] [OAc] (black square), [BMIM+] [OAc] (red circle), and [OMIM+] [OAc] (blue triangle) for (a) Orientational tetrahedral order (q) and (b) Translational tetrahedral order (Sk).

Close modal
FIG. 12.

Probability density distribution of the Orientational tetrahedral order parameter (q) for [OMIM+] [OAc] for 25% (black solid line), 50% (red dashed line), and 75% (blue dotted line) IL concentration and for bulk water (green dashed-dotted line).

FIG. 12.

Probability density distribution of the Orientational tetrahedral order parameter (q) for [OMIM+] [OAc] for 25% (black solid line), 50% (red dashed line), and 75% (blue dotted line) IL concentration and for bulk water (green dashed-dotted line).

Close modal

To analyze the dynamical properties of the IL-water binary mixtures, the self-diffusion coefficients of the cation, anion, and water were calculated from the mean squared displacement (MSD) obtained from the NVE simulations. The resulting diffusion coefficients (Figure 13) are compared to previous simulation work in the literature indicating excellent agreement for comparable systems.27 It is seen that the cation, anion, and water molecules exhibit an exponential decrease in diffusivities with an increase in IL concentration. Furthermore, we can demarcate the decrease in diffusivities into two clear regimes, with a transition occurring at around 75% IL concentration. This demarcation region is observed for all three ILs studied, with the first regime characterized by a drastic decrease in diffusivities as the IL concentration is increased from 25% to 75%. The second region, which corresponds to a further increase in IL concentration from 75% to 100%, is characterized by a gradual decrease in the diffusivities. Similar to the trend in densities, the trend in diffusivities can also be explained based on the changing local structural order of the ILs. Upon initial addition of water to a pure IL, there is a higher ordering in the IL, with a strong polar network present due to the strong electrostatic interactions between the anions and water. There is also a non-polar network due to the tail aggregation, whose effect increases with an increase in cation tail length. Both these effects lead to an ordered structure in the IL s, resulting in low diffusivity values. As the IL concentration is decreased below 75%, a transition occurs, the IL network is disrupted as the water concentration increases and the penetrating water disrupts the ordered IL network. The subsequent addition of water leads to the exponential rise in diffusivities as the system approaches a pure aqueous environment.

FIG. 13.

Effect of dilution on diffusivities for cations (a), anions (b), and water (c). Our results are shown in solid symbols for [EMIM+] [OAc] (square), [BMIM+] [OAc] (circle), and [OMIM+] [OAc] (triangle). Previous simulation work25 is shown in open symbols for [EMIM+] [OAc] (open square) and [BMIM+] [OAc] (open circle).

FIG. 13.

Effect of dilution on diffusivities for cations (a), anions (b), and water (c). Our results are shown in solid symbols for [EMIM+] [OAc] (square), [BMIM+] [OAc] (circle), and [OMIM+] [OAc] (triangle). Previous simulation work25 is shown in open symbols for [EMIM+] [OAc] (open square) and [BMIM+] [OAc] (open circle).

Close modal

A clue to the observed impact of water on the IL is seen when evaluating the impact of the alkyl tail length on the cation diffusivity as the amount of water is modified. At 100% IL, the cation tail length does not have a significant impact on diffusivities, while at 25% IL, the cation diffusivity decreases with increasing tail length, with the transition occurring at 75%, similar to the transitions in the densities. As the concentration of water increases in the system and the IL concentration drops below 75%, the water molecules begin to disrupt the IL network. However, owing to the larger tail lengths, there is an enhanced ordering in the larger alkyl tail ILs driven by a greater number of hydrophobic interactions, which is observed even at very low IL concentrations. The reduced diffusion coefficients for dilute ILs with larger hydrophobic tails are hypothesized to result in slower overall system dynamics with respect to ILs containing smaller hydrophobic tails.

Another interesting observation from the diffusivity results is that for 100% IL, the smaller anion diffuses slower than the much larger cation. This phenomenon has been reported previously in experimentally determined diffusion coefficients30,65,66 and has been attributed to the strong electrostatic interaction between the anions and cations in these systems, where the ions in the systems do not diffuse as single ions but with some amount of ion pairing, hence rendering the effective size of the anions larger than its actual size. The ratio of the anion diffusivity to cation diffusivity for [EMIM+] [OAc] is seen to be ∼0.9, which agrees with experimental observations for this IL.65,67 This ratio decreases with an increase in cation tail length resulting in a ratio of anion diffusivity to cation diffusivity of 0.8 for [BMIM+] [OAc] and to 0.6 for [OMIM+] [OAc]. It is hypothesized that as the length of the cation tail increases, the strength of the cation-anion interactions increases, hence enhancing the pairing of the anion with the cation, which leads to a larger difference between the anion and cation diffusion coefficients. This increase in cation-anion interaction can be seen in Figure 4 where the strength of the interactions increases with an increase in cation tail length, which is reflected in the increasing amplitude as the cation tail length increases. This observation underlies the importance of the alkyl tail length on anion diffusivities.

To further probe the dynamical properties of IL-water binary mixtures, we performed calculations to determine the viscosity of the various solutions. The viscosity was calculated using the Stokes-Einstein equation,68 and despite the highly approximate nature of the calculations (e.g., spherical symmetry, assuming the radius of gyration is approximate to the hydrodynamic radius) they show decent agreement with experimental results69 (Figure 5 in the supplementary material). From the figure, it is observed that the viscosity follows the same trends as that seen in the density and diffusion coefficients, which is not surprising given the use of the diffusion coefficient in the calculation of the viscosity. The viscosity is observed to gradually increase with IL concentration until a concentration of 75%. Beyond 75%, due to the higher ordering of the IL network, the viscosity increases exponentially, which also contributes to the observed slower dynamics of the system.

In this study, we performed a systematic investigation to compute the structural and dynamical properties of three imidazolium based ionic liquids at various water dilutions using fully atomistic molecular dynamics simulations in order to elucidate the effect of water dilution and cation alkyl tail length on the binary mixture characteristics. The results indicate a non-linear transition in system behavior, wherein there is strong ordering in the local structure of IL at high IL concentrations and a breakdown of the structure at low IL concentrations due to the intrusion of water molecules, with a transition occurring at around 75% IL w/w. The computed results also showed good agreement with the literature.

This behavior was demonstrated through the computed density values, which exhibited a maximal density at around 75% IL water content in an IL-water binary mixture. The self-diffusion coefficients were also computed for the cation, anion, and water molecules, and for all the ILs studied, our values showed excellent agreement with the literature. Two distinct regimes were observed for the diffusion values, wherein at low concentrations, they decreased exponentially with increasing IL concentration until a concentration of 75% IL, followed by a more gradual decrease. Viscosity calculations revealed similar trends to the diffusivities.

Through further analysis of the local structure of the ILs in water using RDFs, clustering analysis, order parameters, and interaction energies between the components of the system, the reason for the non-linear behavior of the IL systems was determined to be the local ordering of the ILs in water. When water is added to pure ILs, there exist strong electrostatic interactions between the anions and water in the system. This causes a structured polar network in the ILs, leading to a rise in the density, with a peak reached at around 75% IL concentration. When more water is added to the system, there is a breakdown in the structure of the ILs, manifested in a decrease in the density of the system. As the cation tail length is increased, the cation tail aggregation increases, causing a structuring in the system due to this non-polar network. Hence, for the longer tails of [OMIM+] [OAc], the density turnover occurs at lower concentrations, and as the tail length increases, the effect of tail aggregation on local ordering becomes more significant. Hence, in IL-water systems, there are two factors that contribute to local ordering in the IL structure, anion-water interactions leading to structuring of the polar domain and cation tail aggregation leading to structuring of the non-polar domain, both of which act in unison to resist the breakdown of the IL structure with the addition of water, resulting in a peak in local ordering at around 75% IL concentration. This results in the non-linear trends seen in the properties of the IL water binary mixtures such as densities, diffusion coefficients, and viscosities.

These results provide significant insight into the effect of cationic tail length and water dilution on ILs, hence taking us a step further towards the ultimate goal of being able to design ILs for cellulose pretreatment in order to efficiently utilize lignocellulosic biomass as a viable source of renewable energy.

See supplementary material for additional information as mentioned in the manuscript.

The authors would like to acknowledge the National Science Foundation (NSF Grant No. CBET-1337044) for making this research possible and the Colorado School of Mines High Performance Computing Center.

The authors declare no competing financial interest.

1.
P. C.
Marr
and
A. C.
Marr
,
Green Chem.
18
,
105
128
(
2016
).
2.
A.
Brandt
,
J.
Grasvik
,
J. P.
Hallett
, and
T.
Welton
,
Green Chem.
15
,
550
583
(
2013
).
3.
T.
Welton
,
Coord. Chem. Rev.
248
,
2459
2477
(
2004
).
4.
G. G.
Eshetu
,
M.
Armand
,
H.
Ohno
,
B.
Scrosati
, and
S.
Passerini
,
Energy Environ. Sci.
9
,
49
61
(
2016
).
5.
M.
Doyle
,
S. K.
Choi
, and
G.
Proulx
,
J. Electrochem. Soc.
147
,
34
37
(
2000
).
6.
D. R.
Macfarlane
,
M.
Forsyth
,
P. C.
Howlett
,
J. M.
Pringle
,
J.
Sun
,
G.
Annat
,
W.
Neil
, and
E. I.
Izgorodina
,
Acc. Chem. Res.
40
,
1165
1173
(
2007
).
7.
Y.
Roman-Leshkov
,
C. J.
Barrett
,
Z. Y.
Liu
, and
J. A.
Dumesic
,
Nature
447
,
982
985
(
2007
).
8.
A.
Boisen
,
T. B.
Christensen
,
W.
Fu
,
Y. Y.
Gorbanev
,
T. S.
Hansen
,
J. S.
Jensen
,
S. K.
Klitgaard
,
S.
Pedersen
,
A.
Riisager
,
T.
Stahlberg
, and
J. M.
Woodley
,
Chem. Eng. Res. Des.
87
,
1318
1327
(
2009
).
9.
A.
Milbrandt
,
A Geographic Perspective on the Current Biomass Resource Availability in the United States
, Tech. Rep NREL/TP-560–39181 (National Renewable Energy Laboratory, 2005).
10.
B.
Yang
and
C. E.
Wyman
,
Biofuels, Bioprod. Biorefin.
2
,
26
40
(
2008
).
11.
D.-S.
Lee
,
S. G.
Wi
,
S. J.
Lee
,
Y.-G.
Lee
,
Y.-S.
Kim
, and
H.-J.
Bae
,
Bioresour. Technol.
158
,
239
247
(
2014
).
12.
N.
Muhammad
,
Y. A.
Elsheikh
,
M. I. A.
Mutalib
,
A. A.
Bazmi
,
R. A.
Khan
,
H.
Khan
,
S.
Rafiq
,
Z.
Man
, and
I.
Khan
,
J. Ind. Eng. Chem.
21
,
1
10
(
2015
).
13.
H.
Kawaguchi
,
T.
Hasunuma
,
C.
Ogino
, and
A.
Kondo
,
Curr. Opin. Biotechnol.
42
,
30
39
(
2016
).
14.
A.
Margeot
,
B.
Hahn-Hagerdal
,
M.
Edlund
,
R.
Slade
, and
F.
Monot
,
Curr. Opin. Biotechnol.
20
,
372
380
(
2009
).
15.
S.
Zhu
,
Y.
Wu
,
Q.
Chen
,
Z.
Yu
,
C.
Wang
,
S.
Jin
,
Y.
Ding
, and
G.
Wu
,
Green Chem.
8
,
325
327
(
2006
).
16.
H.
Wang
,
G.
Gurau
, and
R. D.
Rogers
,
Chem. Soc. Rev.
41
,
1519
1537
(
2012
).
17.
C.
Li
,
B.
Knierim
,
C.
Manisseri
,
R.
Arora
,
H. V.
Scheller
,
M.
Auer
,
K. P.
Vogel
,
B. A.
Simmons
, and
S.
Singh
,
Bioresour. Technol.
101
,
4900
4906
(
2010
).
18.
V. B.
Agbor
,
N.
Cicek
,
R.
Sparling
,
A.
Berlin
, and
D. B.
Levin
,
Biotechnol. Adv.
29
,
675
685
(
2011
).
19.
S.
Vasheghani Farahani
,
Y.-W.
Kim
, and
C. A.
Schall
,
Catal. Today
269
,
2
8
(
2016
).
20.
A. M.
da Costa Lopes
and
R.
Bogel-Łukasik
,
ChemSusChem
8
,
947
965
(
2015
).
21.
N. M.
Konda
,
J.
Shi
,
S.
Singh
,
H. W.
Blanch
,
B. A.
Simmons
, and
D.
Klein-Marcuschamer
,
Biotechnol. Biofuels
7
,
1
11
(
2014
).
22.
D. A.
Case
,
T. A.
Durden
,
T. E.
Cheatham
III
,
C. L.
Simmerling
,
J.
Wang
,
R. E.
Duke
,
R.
Luo
,
R. C.
Walker
,
W.
Zhang
,
K. M.
Merz
,
B.
Roberts
,
S.
Hayik
,
A.
Roitberg
,
G.
Seabra
,
J.
Swails
,
A. W.
Goetz
,
I.
Kolossvary
,
K. F.
Wong
,
F.
Paesani
,
J.
Vanicek
,
R. M.
Wolf
,
J.
Liu
,
X.
Wu
,
S. R.
Brozell
,
T.
Steinbrecher
,
H.
Gohlke
,
Q.
Cai
,
X.
Ye
,
J.
Wang
,
M.-J.
Hsieh
,
G.
Cui
,
D. R.
Roe
,
D. H.
Mathews
,
M. G.
Seetin
,
R.
Salomon-Ferrer
,
C.
Sagui
,
V.
Babin
,
T.
Luchko
,
S.
Gusarov
,
A.
Kovalenko
, and
P. A.
Kollman
,
AMBER 12
(
University of California
,
San Francisco
,
2012
).
23.
V. S.
Bharadwaj
,
T. C.
Schutt
,
T. C.
Ashurst
, and
C. M.
Maupin
,
Phys. Chem. Chem. Phys.
17
,
10668
10678
(
2015
).
24.
T.
Schutt
,
V.
Bharadwaj
,
G. A.
Hegde
,
A.
Johns
, and
C. M.
Maupin
, “
In silico insights into the solvation characteristics of the ionic liquid 1-methyl-triethoxy-3-ethylimidazolium acetate for cellulosic biomass
,”
Phys. Chem. Chem. Phys.
(to be published).
25.
T. C.
Schutt
,
G. A.
Hegde
,
V. S.
Bharadwaj
,
A. J.
Johns
, and
C. M.
Maupin
,
Impact of Water-dilution on the Biomass Solvation Properties of an Oligo (Ethoxy) Functionalized Imidazolium Acetate Ionic Liquid
(unpublished).
26.
J. B.
Binder
and
R. T.
Raines
,
Proc. Natl. Acad. Sci. U. S. A.
107
,
4516
4521
(
2010
).
27.
D.
Ghoshdastidar
and
S.
Senapati
,
J. Phys. Chem. B
119
,
10911
10920
(
2015
).
28.
A. A.
Niazi
,
B. D.
Rabideau
, and
A. E.
Ismail
,
J. Phys. Chem. B
117
,
1378
1388
(
2013
).
29.
D.
Ghoshdastidar
and
S.
Senapati
,
Soft matter
12
,
3032
3045
(
2016
).
30.
R. C.
Remsing
,
R. P.
Swatloski
,
R. D.
Rogers
, and
G.
Moyna
,
Chem. Commun.
2006
,
1271
1273
.
31.
R.
Parthasarathi
,
K.
Balamurugan
,
J.
Shi
,
V.
Subramanian
,
B. A.
Simmons
, and
S.
Singh
,
J. Phys. Chem. B
119
,
14339
14349
(
2015
).
32.
Y.
Chen
,
Y.
Cao
,
Y.
Zhang
, and
T.
Mu
,
J. Mol. Struct.
1058
,
244
251
(
2014
).
33.
S.
Stevanovic
,
A.
Podgorsek
,
A. A. H.
Padua
, and
M. F. C.
Gomes
,
J. Phys. Chem. B
116
,
14416
14425
(
2012
).
34.
E.
Quijada-Maldonado
,
S.
van der Boogaart
,
J. H.
Lijbers
,
G. W.
Meindersma
, and
A. B.
de Haan
,
J. Chem. Thermodyn.
51
,
51
58
(
2012
).
35.
C. A.
Hall
,
K. A.
Le
,
C.
Rudaz
,
A.
Radhi
,
C. S.
Lovell
,
R. A.
Damion
,
T.
Budtova
, and
M. E.
Ries
,
J. Phys. Chem. B
116
,
12810
12818
(
2012
).
36.
W.
Guan
,
X.-X.
Ma
,
L.
Li
,
J.
Tong
,
D.-W.
Fang
, and
J.-Z.
Yang
,
J. Phys. Chem. B
115
,
12915
12920
(
2011
).
37.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
,
J. Chem. Phys.
79
,
926
935
(
1983
).
38.
S. V.
Sambasivarao
and
O.
Acevedo
,
J. Chem. Theory Comput.
5
,
1038
1050
(
2009
).
39.
J. M.
Wang
,
R. M.
Wolf
,
J. W.
Caldwell
,
P. A.
Kollman
, and
D. A.
Case
,
J. Comput. Chem.
25
,
1157
1174
(
2004
).
40.
L.
Martinez
,
R.
Andrade
,
E. G.
Birgin
, and
J. M.
Martinez
,
J. Comput. Chem.
30
,
2157
2164
(
2009
).
41.
J. G.
Slingsby
,
S.
Vyas
, and
C. M.
Maupin
,
Mol. Simul.
41
,
1449
1458
(
2015
).
42.
C. L.
Kinsinger
,
Y.
Liu
,
F. L.
Liu
,
Y.
Yang
,
S.
Seifert
,
D. M.
Knauss
,
A. M.
Herring
, and
C. M.
Maupin
,
J. Phys. Chem. C
119
,
24724
24732
(
2015
).
43.
Y.
Liu
,
S. V.
Sambasivarao
,
J. L.
Horan
,
Y.
Yang
,
C. M.
Maupin
, and
A. M.
Herring
,
J. Phys. Chem. C
118
,
854
863
(
2014
).
44.
S. V.
Sambasivarao
,
Y.
Liu
,
J. L.
Horan
,
S.
Seifert
,
A. M.
Herring
, and
C. M.
Maupin
,
J. Phys. Chem. C
118
,
20193
20202
(
2014
).
45.
J.-P.
Ryckaert
,
G.
Ciccotti
, and
H. J.
Berendsen
,
J. Comput. Phys.
23
,
327
341
(
1977
).
46.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
,
J. Mol. Graphics Modell.
14
,
33
38
(
1996
).
47.
J.
Schindelin
,
I.
Arganda-Carreras
,
E.
Frise
,
V.
Kaynig
,
M.
Longair
,
T.
Pietzsch
,
S.
Preibisch
,
C.
Rueden
,
S.
Saalfeld
,
B.
Schmid
,
J. Y.
Tinevez
,
D. J.
White
,
V.
Hartenstein
,
K.
Eliceiri
,
P.
Tomancak
, and
A.
Cardona
,
Nat. Methods
9
,
676
682
(
2012
).
48.
Y.
Wang
and
G. A.
Voth
,
J. Phys. Chem. B
110
,
18601
18608
(
2006
).
49.
Y.
Ji
,
R.
Shi
,
Y.
Wang
, and
G.
Saielli
,
J. Phys. Chem. B
117
,
1104
1109
(
2013
).
50.
N.
Galamba
,
J. Phys. Chem. B
117
,
2153
2159
(
2013
).
51.
N.
Giovambattista
,
P. G.
Debenedetti
,
F.
Sciortino
, and
H. E.
Stanley
,
Phys. Rev. E
71
,
061505
(
2005
).
52.
Z.
Yan
,
S. V.
Buldyrev
,
P.
Kumar
,
N.
Giovambattista
,
P. G.
Debenedetti
, and
H. E.
Stanley
,
Phys. Rev. E
76
,
051201
(
2007
).
53.
D. R.
Nutt
and
J. C.
Smith
,
J. Am. Chem. Soc.
130
,
13066
13073
(
2008
).
54.
D.
Nayar
and
C.
Chakravarty
,
Phys. Chem. Chem. Phys.
15
,
14162
14177
(
2013
).
55.
F.
Sedlmeier
,
D.
Horinek
, and
R. R.
Netz
,
J. Am. Chem. Soc.
133
,
1391
1398
(
2011
).
56.
E.
Duboué-Dijon
and
D.
Laage
,
J. Phys. Chem. B
119
,
8406
8418
(
2015
).
57.
J. R.
Errington
and
P. G.
Debenedetti
,
Nature
409
,
318
321
(
2001
).
58.
P.-L.
Chau
and
A.
Hardwick
,
Mol. Phys.
93
,
511
518
(
1998
).
59.
C.
Caleman
and
D.
van der Spoel
,
J. Chem. Phys.
125
,
154508
(
2006
).
60.
P.
Jedlovszky
,
M.
Mezei
, and
R.
Vallauri
,
Chem. Phys. Lett.
318
,
155
160
(
2000
).
61.
P.
Jedlovszky
,
J. Chem. Phys.
111
,
5975
5985
(
1999
).
62.
W.
Jiang
,
Y.
Wang
, and
G. A.
Voth
,
J. Phys. Chem. B
111
,
4812
4818
(
2007
).
63.
T.
Yan
,
Y.
Wang
, and
C.
Knox
,
J. Phys. Chem. B
114
,
6905
6921
(
2010
).
64.
T.
Yan
,
Y.
Wang
, and
C.
Knox
,
J. Phys. Chem. B
114
,
6886
6904
(
2010
).
65.
D. T.
Bowron
,
C.
D’Agostino
,
L. F.
Gladden
,
C.
Hardacre
,
J. D.
Holbrey
,
M. C.
Lagunas
,
J.
McGregor
,
M. D.
Mantle
,
C. L.
Mullan
, and
T. G. A.
Youngs
,
J. Phys. Chem. B
114
,
7760
7768
(
2010
).
66.
A.
Menjoge
,
J.
Dixon
,
J. F.
Brennecke
,
E. J.
Maginn
, and
S.
Vasenkov
,
J. Phys. Chem. B
113
,
6353
6359
(
2009
).
67.
M. E.
Ries
,
A.
Radhi
,
A. S.
Keating
,
O.
Parker
, and
T.
Budtova
,
Biomacromolecules
15
,
609
617
(
2014
).
68.
V. S.
Bharadwaj
,
N. M.
Eagan
,
N. M.
Wang
,
M. W.
Liberatore
, and
C. M.
Maupin
,
ChemPhysChem
16
,
2810
2817
(
2015
).
69.
S.
Fendt
,
S.
Padmanabhan
,
H. W.
Blanch
, and
J. M.
Prausnitz
,
J. Chem. Eng. Data
56
,
31
34
(
2011
).

Supplementary Material