Molecular insight into the microstructure and microscopic dynamics of pyridinium ionic liquids with different alkyl chains based on temperature response

Abstract

The temperature dependence of the microstructure and microscopic dynamics of three N-alkylpyridinium tetrafluoroborate ionic liquids, namely, [BPy][BF₄], [HPy][BF₄], and [OPy][BF₄], is investigated via molecular dynamics simulation. Under a temperature change, the microstructure changes of [OPy][BF₄] with a longer alkyl chain show a slight difference from that of [BPy][BF₄] and [HPy][BF₄] in terms of the radial and spatial distribution functions. A temperature increase is found to be advantageous to the aggregation of the polar regions as well as the nonpolar regions in the pyridinium ionic liquids. However, the properties of the microscopic dynamics of the three ionic liquids regularly change under conditions with temperature change. The results show a decreased probability of ions to reverse the direction of their translation motion because the negative parts of the velocity autocorrelation functions increase with the temperature. The temperature dependence of rotational dynamics is more obvious than that of translation. A rise in temperature causes rotational autocorrelation functions to decay quickly, and it considerably reduces relaxation time. With an increase in alkyl chain length, the rotational correlation improves and the relaxation time decreases. The time correlation functions show that the association dynamics of nonpolar regions is faster than that of polar regions over time. The association time is short for long-tailed ionic liquids.

---

1. Introduction

Ionic liquids (ILs) are compounds composed of cations and anions, which are liquid at room temperature.^1,2 ILs are characterized by low vapor pressure, low volatility, low melting point, wide liquid range, wide electrochemical window, and so on. As a new type of green solvent, ILs are receiving increasing research attention.^3–5 According to the articles published in recent years, the research on ILs mainly focuses on organic reactions, materials chemistry, polymer chemistry, analytical chemistry, and the chemical industry.^6,7 ILs are also applied to other fields, such as energy, environment, bioscience, photoelectric information, and aerospace materials.^8,9 Owing to the many applications of ILs in different fields, their chemical and physical properties merit further attention.^10,11 Different structures are known to have different properties.¹² The study on the relationship between the microstructure and macroscopic properties of ILs has achieved significant progress; however, some of the data and results are conflicting.^13,14 Thus, we still lack in-depth knowledge of the microscopic structure of ILs and the interaction mechanism of cations and anions.^15,16 Consequently, the research into the microscopic structure and properties of ILs, especially the effects of their dynamics properties on solute–solvent interaction, is a challenging undertaking.^12,17

With the deepening of the research into ILs, we should establish the structure and properties of ILs at the molecular level. The microstructure and thermodynamic properties of ILs can be analyzed and predicted via molecular dynamics (MD) simulation.^12,18–22 A large body of research has explored the microstructure of imidazolium ILs.^21,22 For pyridinium ILs, studies have been dominantly conducted on the experimental aspect.^23–25 For example, thermodynamic properties including densities and viscosities,²⁵ have been measured at certain temperatures. However, a limited amount of research has studied the structures and properties of pyridinium ILs via MD simulation in different temperatures.²⁶ In the present work, MD simulation was carried out to study the pure systems of N-butylpyridinium tetrafluoroborate ([BPy][BF₄]), N-hexylpyridinium tetrafluoroborate ([HPy][BF₄]), and N-octylpyridinium tetrafluoroborate ([OPy][BF₄]). The effects of temperature on the local microstructure and dynamics of the three ILs were compared by analyzing their radial distribution functions, spatial distribution functions, and correlation functions. The relationship between the microstructure and the properties was also explored.

2. Simulation methodology

2.1 Force field parameters

The modified OPLS-AA force field developed by Lopes and Pádua was used for the cations of the ILs,²⁷ as shown in eqn (1).where V_total is the total energy of the system, which includes bond length, bond angle, dihedral angle, van der Waals (VDW) interactions, and electrostatic interactions. The nonbonded interactions comprise Lennard-Jones and coulombic potential terms. The force constants for the bonded interactions are expressed as k_b, k_θ, and k_φ. The bond length, bond angle, and dihedral angle for the equilibrium structure are expressed as r_eq, θ_eq, and δ, respectively. q_i is the charge of atom i, and r_ij is the distance between atoms i and j. ε_ij and δ_ij are the LJ parameters for different atoms. The force field parameters of anion [BF₄]⁻ are taken from the literature.²⁸ Fig. 1 shows the labels of the atoms for the pyridinium cations and anions.

Fig. 1 Molecular structure of ILs in this work: [BPy][BF₄] (red ellipse), [HPy][BF₄] (green ellipse), and [OPy][BF₄] (blue ellipse).

2.2 Simulation details

MD simulations were run using the GROMACS 4.5.6 program.^29,30 Starting configurations containing 256 ILs were subjected to energy minimization, followed by a 10 ns equilibration in the isothermal–isobaric ensemble. On the basis of the results, each system was equilibrated for 3 ns at a setting temperature that was changed from 293 K to 373 K at an interval of 10 K. The last 2 ns trajectories were taken to calculate the density. During this period, the temperature was maintained via a V-rescale thermostat³¹ while the pressure was held at 1 bar with a Berendsen barostat.³² A configuration of the system at the average density was taken, and simulations were then run for 4 ns in the canonical ensemble. Without pressure coupling, the V-rescale method was used to maintain the temperature of the system. The particle mesh Ewald method was applied to the long electrostatic interactions, and the bond length was constrained with the LINCS algorithm.^33,34 The integration time step was 2 fs. The cutoff radius for VDW and the Coulomb interaction were set to 1.5 nm. The standard periodic boundary condition was used in all simulations. The structure analysis was performed using the last 1 ns trajectory. The rotation and association dynamics were sampled over 2 ns of production runs.

3. Results and discussion

3.1 Evaluation of force fields

To verify the accuracy of the force field parameters,³⁵ we obtained the densities of [BPy][BF₄], [HPy][BF₄], and [OPy][BF₄] under different temperatures. In Fig. 2, the hollow part represents the simulated values in the present work, and the solid part represents the experimental values from the literature.^23–25 The densities of the system decrease with increase in temperature and present an almost linear relationship. With the increase of the cationic alkyl chain length, the density values decrease because the distance between particles and the macroscopic volume change. Moreover, the gap of the simulation and experiment values is narrowed with the increase of the chain length. The simulation results are in good agreement with the experimental data, with the percentage error being smaller than 2.0%. The chosen force fields of the ILs are thus accurate enough to use. In addition, the ILs as green solvents can be designed with suitable densities that are compatible with those in other experiments by changing their structure appropriately.³⁶

Fig. 2 Densities of three pyridinium ILs as the functions of temperature; the experimental values (solid) are taken from ref. 23–25.

3.2 Microstructure

As described in the work of Voth³⁷ and our previous study,³⁸ the N-alkylpyridinium tetrafluoroborate ILs can be divided into two regions: the “head” or the polar region, which is composed of a pyridine ring of cation and anion [BF₄]⁻, and the “tail” part or the nonpolar region, which is composed of the alkyl chain of cations. If this description of the heterogeneous structure of N-alkylpyridinium ILs is correct, the different interactions might cause variations in the effects of the temperature dependence of the polar and nonpolar regions.

3.2.1 Radial distribution functions. The radial distribution functions (RDFs) describe the characteristics of short-range order and can reflect the microscopic structure of a material.³⁹ Generally, RDFs are used to describe the structure of ILs effectively.⁴⁰ In this section, we discuss the influence of temperature on the local microstructure of pyridinium ILs based on RDFs. To clearly recognize the temperature dependence effect, we compare in detail the changes in the site–site RDFs of the CTn–CTn of tails (n = 4, 6, 8, corresponding to the ending carbon atom in the alkyl chain of the three ILs) and HA1–B of heads at different temperatures.

As shown in Fig. 3, the RDFs between the terminal carbon atoms of the alkyl chain were calculated at different temperatures. With the increase of the temperature, the RDFs show that the peaks of CT4–CT4 of [BPy][BF₄] and CT6–CT6 of [HPy][BF₄] are reduced. This result indicates that the aggregation degree of the ending carbon atoms decreases with the weakening of their interactions. The peak of [HPy][BF₄] is higher than the peak of [BPy][BF₄]. The peaks of CT8–CT8 from [OPy][BF₄] rise progressively with the increase in the temperature, and the peaks are relatively wide. This indicates that the CT8–CT8 aggregation degree is enhanced. Moreover, all the positions of the peaks show almost no changes at r = 0.42 nm. According to the above results, the change in the nonpolar region of [OPy][BF₄] is obviously different from those of [BPy][BF₄] and [HPy][BF₄] with the increase in temperature.

Fig. 3 Site–site RDFs for the terminal carbon atom (CT4, CT6 and CT8) of the alkyl chains for [BPy][BF₄], [HPy][BF₄] and [OPy][BF₄] respectively.

Fig. 4 shows the site–site RDFs of the F atoms of BF₄ around the HA1 atoms of the cation. The position of the first peak is at the maximum at r = 0.25 nm, which is the characteristic hydrogen bond distance of the C–H⋯F pairs.⁴¹ This result indicates a strong interaction between HA1 and F at this position. The position of the second peak is at 0.44 nm. With the increase in the temperature, the first peaks of the RDFs decrease. The changes of the second peaks are not obvious. With the increase of the alkyl chain length, the peaks increase in magnitude. According to Fig. 4, all the changes in the interaction between HA1 and F are similar when the temperature rises; the same is not true for CTn. With the increase in the temperature, the intensity of the first peaks decreases, possibly because of the temperature dependence of the hydrogen bond.⁴²

Fig. 4 Site–site RDFs of the F atoms of BF₄ around the HA1 atoms of the cation.

To comprehensively study the temperature dependence of the microstructure of the polar region, we analyze the RDFs of the center-of-mass (COM) between the pyridine ring of the cation and the anion.⁴³ As shown in Fig. 5, the COM RDFs are observed at 293 K and 373 K. More COM RDFs are presented in Fig. S2 of the ESI.† The position of the highest peak is at 0.52 nm, which is the first shell. This position indicates that the maximum probability of the anion around the head is at approximately 0.52 nm. The second peak at 1.05 nm, which is the second shell, is relatively weak. With the increase in the temperature, the peaks of [BPy][BF₄] and [HPy][BF₄] are reduced, and the second peaks basically remain unchanged. This result suggests that the rise of the temperature weakens the interaction between the cations and the anions of these two ILs. By contrast, the maximum peak of [OPy][BF₄] increases, and the second peak slightly increases. The structural change of [OPy][BF₄] obviously shows a slight difference from those of [BPy][BF₄] and [HPy][BF₄] with the change in the temperature. This difference is explained in the next section detailing the coordination number (CN).

Fig. 5 Center-of-mass RDFs for anions around the pyridine ring of cations.

Fig. 5 shows that with the increase of the alkyl chain length, the height of the first peak of the RDFs increases, and the width widens. This condition illustrates that the interaction between the cation and the anion strengthened because of the increase in the alkyl chain length. However, all the positions of the peaks show almost no changes, where the atomic number density being significantly larger than the average.

3.2.2 Coordination number. The CN is the atomic number with the central atoms coordinating.⁴⁴ CNs can help us understand the microstructures of the pyridinium ILs. Generally, the CNs of the first shell in different systems are gained by integrating the RDFs into the location of the first minimum.⁴⁵ The calculation formula of the CNs can be expressed as eqn (2).where ρ is the number density and r_shell is the first minimum in the RDFs.

Fig. 6(a) shows the CNs of the terminal carbon atom CTn of the tail for the three ILs. With the increase in the temperature, the CNs of CT4–CT4 and CT6–CT6 decrease, and the CNs of [HPy][BF₄] become larger than those of [BPy][BF₄] at the same temperature. However, the CNs of CT8–CT8 increase obviously. As shown in Fig. 6(b), with the increase in the temperature, the CNs of the F atoms around the HA1 atom decrease, and the curve of [BPy][BF₄] falls rapidly. According to Fig. 6(a) and (b), with the change in the temperature, the CNs for HA1–F and CT8–CT8 of [OPy][BF₄] obviously differ from those of [BPy][BF₄] and [HPy][BF₄]. This result indicates that the longer tails of [OPy][BF₄] may induce changes in the local structure because the three ILs show a similar composition, except for the different alkyl chains.

Fig. 6 Coordination numbers of tails CTn–CTn (a), HA1–F (b), and pyridine ring–anion (c) for the three ILs under different temperatures.

To further realize this change induced by the temperature, we calculate the CNs for the pyridine ring–anion (heads) on the basis of the RDFs of COM of the pyridine ring–anion. As shown in Fig. 6(c), the coordination ability of the head of [BPy][BF₄] is the strongest in all temperature ranges. With the increase in the temperature, the CNs of the heads of [BPy][BF₄] and [HPy][BF₄] decrease, whereas the CNs of [OPy][BF₄] increase. That is, the CNs of both heads and tails of [OPy][BF₄] increase with the rise in the temperature, while the volume expands and the density decreases.

Generally, as the alkyl chain lengths increase, the number density should decrease correspondingly to reduce coordination ability gradually. The abnormality of [OPy][BF₄] is explained as follows. [OPy][BF₄] has a relatively long alkyl chain, which favors the formation of large clusters.⁴⁶ For [OPy][BF₄], the rise in temperature causes the aggregation degree of the polar and nonpolar regions to increase, which in turn results in the increase of the RDFs and CNs. Meanwhile, the gaps between the large clusters increase with the rise in the temperature. Therefore, the volume of the system still expands, and the density decreases.

3.2.3 Spatial distribution functions. Spatial distribution functions (SDFs) can intuitively present the three-dimensional distribution of each system,⁴⁷ which denotes the distribution probability of the particles around the center particle b. We derive the SDFs of the anions around the cation and the CTn around the cation of the three systems under different temperatures. As shown in Fig. 7, the SDFs of the three ILs are drawn using the VMD software.⁴⁸ Given that the spatial density in the figures is relative, we adjust the isovalue to its improper value so that the figure of the SDFs could be clearly displayed.

Fig. 7 Spatial distributions of anions around cations in red color and those of CTn around cations in yellow color at 293 K (left, the isovalues are 35 for red and 45 for yellow) and 373 K (right, the isovalues are 20 for red and 25 for yellow).

The RDFs and SDFs show that the anions are mainly distributed around the cationic pyridine ring. This distribution is due to the strong electrostatic interactions and hydrogen bond between HA1 on the rings and F on [BF₄]⁻. CTn (n = 4, 6, 8) is mainly distributed around the tail of the alkyl chain for the VDW interaction. Under the same low temperature, the anions are distributed widely with the increase of the alkyl chain length. Under the same high temperature, the aggregation of the polar and nonpolar regions could be observed clearly with the increase of the chain length.

Furthermore, the probabilities of the anions around the pyridine ring increase with the rise in the temperature, as observed by comparing the SDFs at 293 K and 373 K in Fig. 7. A similar phenomenon could be observed for the CTn around the tail. This temperature effects are not only observed for the long chain [OPy][BF₄] but also for the other ILs. More SDFs of the anions around the pyridine ring are shown in Fig. S3 of ESI.† These changes cannot be observed from the site–site RDFs; thus, SDFs can present the overall structural changes of ILs. In addition, the rise in the temperature causes the cations and anions to translate and rotate rapidly, which in turn reduces viscosity.^24,49 Specifically, a high temperature within the ranges of this study is advantageous to the aggregation of polar regions as well as nonpolar regions, as indicated by the following microscopic dynamics results. This unusual characteristic has not been observed in conventional solvents, which could help us understand the particular solute–solvent behavior in ILs.⁵⁰ The aggregation behavior with increasing temperature maybe results in different local solvation structure in polar and nonpolar regions owing to different solute–solvent interactions. In ILs, the appropriate temperature will be important to the diffusion, dissolution and reaction of solute.

3.3 Microscopic dynamics

3.3.1 Velocity autocorrelation functions. Velocity autocorrelation functions (VACFs) can characterize molecular diffusion properties and describe the velocity correlation degree of the center mass.⁵¹ As shown in Fig. 8, the VACFs of the cations and anions of [HPy][BF₄] vary with time under different temperatures. The VACFs drop rapidly, decay slightly, and then increase to zero gradually. This trend shows that the correlation of velocity to the initial moment changes over time. Noticeably, the motions present a negative part. As depicted in Fig. 8, the increase of the temperature causes the diminishing of the negative part of the VACFs to a certain degree. This condition indicates the decreased probability for the particle to reverse the direction of its motion.⁵² The temperature responses for the translation movements of [BPy][BF₄] and [OPy][BF₄] are similar to that of [HPy][BF₄] (see Fig. S4 in ESI†). The VACFs change quickly in the order of one picosecond. For the cations, no correlation is observed after 1.3 ps, which is the relaxation time since the VACFs fluctuate to zero. The relaxation time of the anions is about 0.8 ps.

Fig. 8 Velocity autocorrelation functions of the cations and anions of [HPy][BF₄] under different temperatures.

As shown in Fig. 9, the attenuation of the VACFs for [BPy]⁺ is the fastest at each temperature. As the alkyl chain length increases, the negative part increases. This result indicates that reversing the direction of the motion in ILs with long alkyl chains is increasingly difficult because of the high viscosities.²⁴ The anionic VACFs of the ILs change slightly at the same temperature because of their small scale (Fig. S5 in ESI†). The temperature shows certain effects on the translational motion. In the next section, we explore rotational movement.

Fig. 9 Comparison of the velocity autocorrelation functions of cations with different alkyl chain lengths at 293 K (a) and 373 K (b).

3.3.2 Rotational autocorrelation functions. We choose the pyridine ring to study how its rotational dynamics is affected by temperature. The rotational autocorrelation functions (RACFs) are shown in Fig. 10, which reflects the temperature dependence of the rate of reorientation or rotational motion of the pyridine ring. As shown in Fig. 10, the decline of the RACFs evolves over time until they finally reach zero. With the increase in the temperature, the RACFs decay rapidly, thus indicating the decrease in the rotational correlation. At the same temperature, the attenuation becomes slow with the increase of the alkyl chain length, thus showing an improved correlation.

Fig. 10 RACFs of perpendicular line in the pyridine ring under different temperatures for the three ILs.

The relaxation times are fitted from the curves in Fig. 10 to a double exponential decay,⁵³ as shown in eqn (3).

where τ₁ and τ₂ represent the relaxation times of ILs and a is a proportion parameter. All variables are listed in Table 1. With the increase in the temperature, the short relaxation time τ₁ and long relaxation time τ₂ decrease. Both τ₁ and τ₂ increase with the alkyl chain length. Thus, as shown in the relaxation times in Table 1, the temperature exerts a great effect on the local rotational dynamics of pure ILs. The result that temperatures make microscopic dynamics faster is similar as the situation adding the molecular solvent.²⁰

Table 1 Changes of the rotational relaxation time and proportion parameter (a) for the pyridine ring of the three ILs with temperature

ILs	τ1 (ps) (a)				τ2 (ps)
	293 K	323 K	353 K	373 K	293 K	323 K	353 K	373 K
[BPy][BF4]	30.89 (0.267)	20.85 (0.306)	10.41 (0.298)	5.37 (0.270)	957.92	432.79	181.61	92.77
[HPy][BF4]	41.28 (0.210)	34.66 (0.301)	20.05 (0.317)	11.44 (0.306)	1607.90	906.31	460.54	206.54
[OPy][BF4]	45.57 (0.197)	38.34 (0.277)	31.02 (0.354)	24.39 (0.403)	2083.70	1158.50	578.04	355.16

---

3.3.3 Time correlation functions. We further study how temperature influences the association dynamics in polar and nonpolar regions by comparing the time correlation functions (TCFs).⁵⁴ A TCF is defined as , where θ(r,t) represents the Heaviside step function.¹⁷ When the two associated particles are combined, θ(r,t) gets a value of 1, and when the two particles depart, the value becomes zero. N_C is the number of selected particles, and N_A is the number of associated particles that are bound to the selected particles. Here, we choose HA1–B to study the association in the polar region. If the distance is within the minimum of the first coordination shell of their RDFs, HA1 and B are deemed associated. Similarly, the association in the nonpolar region is determined by comparing the distance of CTn–CTn (n = 4, 6, 8) with their RDF minimum at the first coordination shell.

The temperature dependence effects on the association dynamics of [HPy][BF₄] in the polar and nonpolar regions are compared in Fig. 11 in terms of TCFs. Other ILs can be found in Fig. S6–S8 of the ESI.† Evolving over time, the TCFs first drop fast and then decay slowly. The long-time decay of HA1–B is extremely slow under low temperature, thus indicating the strong interactions between the cations and the anions. As the temperature increases, the TCFs of the three pyridinium ILs decay quickly and gradually approach zero within an adequate time. Conversely, the association time of the nonpolar regions is shorter than that of the polar regions. This result indicates a relatively fast association dynamics for the alkyl chain because the tails are associated because of the VDW interactions, which are weaker than the Coulomb interactions in the polar regions.⁵⁵

Fig. 11 Time correlation functions of HA1–B (a) and ending carbon CT6–CT6 (b) for [HPy][BF₄] at different temperatures.

As shown in Fig. 12, with the increase of the alkyl chain length, the tendency of the function to decay to zero accelerates. The attenuation of [OPy][BF₄] is the fastest. The comparison of Fig. 12(a) and (b) shows that the attenuation time is longer at 293 K than at 373 K. The TCFs of the three ILs are not zero in 1 ns at 293 K, whereas all the functions decay to zero in about 0.5 ns at 373 K. Hence, temperature exerts a significant effect on ion association. Although the temperature increase leads to its reduction within the average association time, the aggregation degree of the polar regions increases for the quick interchange of anions and cations among the ion pairs.

Fig. 12 TCFs of HA1–B of three ILs with different alkyl chain lengths at 293 K ((a) the solid symbol) and 373 K ((b) the hollow symbol).

4. Conclusion

A systematic study is conducted on the local structure and the microscopic dynamics of three types of ILs. The influence of temperature on the structures of ILs is obvious. For [OPy][BF₄] with one longer alkyl chain, some particularity is observed in its structural changes when the temperature rises. Such particularity is reflected by the RDFs and SDFs. Most of the site–site RDFs decrease with the increase in the temperature, except for the tail of [OPy][BF₄]. The SDFs indicate an increase in aggregation at a high temperature for both the polar and the nonpolar regions in the pyridinium ILs. The anions are mainly distributed around the pyridine rings, and the tails of the cations gather together.

In microscopic dynamics, as the temperature increases, the molecular thermal motion becomes fast and the correlation declines. The resulting VACFs show that a high temperature and long alkyl chain cause ions to easily reverse the direction of their motion. Over time, the rotation of the pyridine rings decays rapidly with the increase in the temperature. This outcome leads to a decrease in the correlation and relaxation time of rotation. We also study how temperature influences the association in the polar and nonpolar regions on the basis of TCFs. With the increase in the temperature, the association dynamics decrease. At the same temperature, the association time of [OPy][BF₄] is the shortest for its longest chains in the three ILs. Under the same system, the association time for the polar region is longer than that for the nonpolar region. These results will be helpful to understand the particular solute–solvent behaviors in ionic liquids.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (No. 21173002, 61475001, 21103040 and 21673003) and the Science Foundation of the Education Committee of Anhui Province of China (No. KJ2010A145).

References

1. R. D. Rogers and K. R. Seddon, Science, 2003, 302, 792–793.
2. Z. Liu, S. Huang and W. Wang, J. Phys. Chem. B, 2004, 108, 12978–12989.
3. Ł. Marcinkowski, A. Kloskowski, J. Czub, J. Namieśnik and D. Warmińska, J. Chem. Thermodyn., 2015, 88, 36–43.
4. T. Pal and R. Biswas, J. Phys. Chem. B, 2015, 119, 15683–15695.
5. N. Deenadayalu, I. Bahadur and T. Hofman, J. Chem. Eng. Data, 2011, 56, 1682–1686.
6. K. C. Badgujar and B. M. Bhanage, Bioresour. Technol., 2015, 178, 2–18.
7. H. Olivier-Bourbigou, L. Magna and D. Morvan, Appl. Catal., A, 2010, 373, 1–56.
8. D. Wei and A. Ivaska, Anal. Chim. Acta, 2008, 607, 126–135.
9. G. Quijano, A. Couvert and A. Amrane, Bioresour. Technol., 2010, 101, 8923–8930.
10. N. Deenadayalu, I. Bahadur and T. Hofman, J. Chem. Thermodyn., 2010, 42, 726–733.
11. V. Govinda, P. M. Reddy, I. Bahadur, P. Attri, P. Venkatesu and P. Venkateswarlu, Thermochim. Acta, 2013, 556, 75–88.
12. R. Hayes, G. G. Warr and R. Atkin, Chem. Rev., 2015, 115, 6357–6426.
13. M. Liang, A. Kaintz, G. A. Baker and M. Maroncelli, J. Phys. Chem. B, 2012, 116, 1370–1384.
14. C. Chiappe and C. S. Pomelli, Phys. Chem. Chem. Phys., 2013, 15, 412–423.
15. P. Henritzi, A. Bormuth, F. Klameth and M. Vogel, J. Chem. Phys., 2015, 143, 164502.
16. G. W. Meindersma, A. R. Hansmeier and A. B. D. Haan, Ind. Eng. Chem. Res., 2010, 49, 7530–7540.
17. P. K. Sahu, S. K. Das and M. Sarkar, J. Phys. Chem. B, 2014, 118, 1907–1915.
18. Y. Wang, M. R. Shimpi, S. Sarman, O. N. Antzutkin, S. Glavatskih, L. Kloo and A. Laaksonen, J. Phys. Chem. B, 2016, 120, 7446–7455.
19. E. J. Maginn, J. Phys.: Condens. Matter, 2009, 21, 373101.
20. M. Sha, H. Dong, F. Luo, Z. Tang, G. Zhu and G. Wu, J. Phys. Chem. Lett., 2015, 6, 3713–3720.
21. G. E. Logotheti, J. Ramos and I. G. Economou, J. Phys. Chem. B, 2009, 113, 7211–7224.
22. K. R. Ramya, P. Kumar and A. Venkatnathan, J. Phys. Chem. B, 2015, 119, 14800–14806.
23. D. Tomida, S. Kenmochi, K. Qiao, T. Tsukada and C. Yokoyama, Fluid Phase Equilib., 2013, 340, 31–36.
24. B. Mokhtarani, A. Sharifi, H. R. Mortaheb, M. Mirzaei, M. Mafi and F. Sadeghian, J. Chem. Thermodyn., 2009, 41, 323–329.
25. M. Larriba, S. García, P. Navarro, J. García and F. Rodríguez, J. Chem. Eng. Data, 2012, 57, 1318–1325.
26. I. V. Voroshylova and V. V. Chaban, J. Phys. Chem. B, 2014, 118, 10716–10724.
27. J. N. A. C. Lopes and A. A. H. Pádua, J. Phys. Chem. B, 2006, 110, 19586–19592.
28. X. Wu, Z. Liu, S. Huang and W. Wang, Phys. Chem. Chem. Phys., 2005, 7, 2771–2779.
29. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
30. H. J. C. Berendsen, D. van der Spoel and R. van Drunen, Comput. Phys. Commun., 1995, 91, 43–56.
31. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
32. H. J. Berendsen, J. V. Postma, W. F. van Gunsteren, A. R. H. J. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
33. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
34. B. Hess, H. Bekker, H. J. Berendsen and J. G. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
35. D. L. Cheung, S. J. Clark and M. R. Wilson, Chem. Phys. Lett., 2002, 356, 140–146.
36. I. Bahadur and N. Deenadayalu, J. Solution Chem., 2011, 40, 1528–1543.
37. W. Jiang, Y. Wang and G. A. Voth, J. Phys. Chem. B, 2007, 111, 4812–4818.
38. G. Zhu, Y. Wang, L. Zhang, Y. Luo, M. Sha and X. Xu, J. Mol. Liq., 2015, 203, 153–158.
39. V. Lesch, Z. Li, D. Bedrov, O. Borodin and A. Heuer, Phys. Chem. Chem. Phys., 2016, 18, 382–392.
40. C. Knight, C. M. Maupin, S. Izvekov and G. A. Voth, J. Chem. Theory Comput., 2010, 6, 3223–3232.
41. H. Sun, B. Qiao, D. Zhang and C. Liu, J. Phys. Chem. A, 2010, 114, 3990–3996.
42. I. Bandrés, R. Alcalde, C. Lafuente, M. Atilhan and S. Aparicio, J. Phys. Chem. B, 2011, 115, 12499–12513.
43. S. Li, G. Feng, J. L. Bañuelos, G. Rother, P. F. Fulvio, S. Dai and P. T. Cummings, J. Phys. Chem. C, 2013, 117, 18251–18257.
44. Y. Chen, N. L. Doltsinis, R. C. Hider and D. J. Barlow, J. Phys. Chem. Lett., 2012, 3, 2980–2985.
45. I. Morrow Timothy and J. Maginn Edward, J. Phys. Chem. B, 2002, 106, 12807–12813.
46. J. Wang, H. Wang, S. Zhang, S. Zhang and Y. Zhao, J. Phys. Chem. B, 2007, 111, 6181–6188.
47. K. Singer, J. V. L. Singer and A. J. Taylor, Mol. Phys., 1979, 37, 1239–1262.
48. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
49. T. Kavitha, P. Attri, P. Venkatesu, R. R. Devi and T. Hofman, J. Phys. Chem. B, 2012, 116, 4561–4574.
50. T. Chan, H. Li and K. Li, J. Phys. Chem. B, 2015, 119, 15718–15728.
51. D. Roy, S. Liu, B. L. Woods, A. R. Siler, J. T. Fourkas, J. D. Weeks and R. A. Walker, J. Phys. Chem. C, 2013, 117, 27052–27061.
52. S. Bouazizi and S. Nasr, J. Mol. Liq., 2011, 162, 78–83.
53. M. E. Siwko, S. J. Marrink, A. H. de Vries, A. Kozubek, A. J. S. Uiterkamp and A. E. Mark, Biochim. Biophys. Acta, Biomembr., 2007, 1768, 198–206.
54. T. Chang, L. Dang, R. Devanathan and M. Dupuis, J. Phys. Chem. A, 2010, 114, 12764–12774.
55. T. A. Lima, V. H. Paschoal, L. F. Faria, M. C. Ribeiro and C. Giles, J. Chem. Phys., 2016, 144, 224504.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra26866h

---