Ultraconfinement of aqueous electrolytic solutions within hydrophilic nanotubes

Abstract

By means of molecular simulations we shed light on the interplay of surface, confinement and salt effects on the structure and dynamics of water and ions highly confined within a hydrophilic silica nanotube. By a local description of the hydrogen bonding network, the structure and the dynamics we disentangle the confinement and the surface effects. In general, the concentration dependence of the structure and the dynamics is strongly exacerbated under high confinement. Thus, the confinement effect enhances the formation of Na–Cl clusters which strongly affect the dynamics of water. We show that the interlayer which corresponds to a minimum in radial density is linked to an unexpected increase in the hydrogen bonds number while dynamically we highlight a subdiffusive regime together with an increase of the interfacial dipolar relaxation time for water molecules as the NaCl concentration increases. Energetically speaking we demonstrate that the cohesion of the HB network decreases significantly under highly restrictive geometry.

---

1 Introduction

Water is ubiquitous in all sectors of society, from drinking to agriculture and from energy supply to industrial manufacturing. With the increasing demand for water and diminishing freshwater sources as a result of environmental degradation, new technologies for water supply have a crucial role to play in addressing the worlds clean water needs in this century. Desalination is in many regards the most promising approach to long term water supply. Seawater desalination using reverse osmosis (RO) membranes has become a common method for countries with direct access to the sea. While this technology is widely used and has proven to be efficient, it remains, however, relatively costly due to the use of high pressure pumps and to the small permeation rate of RO membranes. In order to improve and better interpret their performances it is compulsory to understand both structural and dynamical properties of confined water and ions into nanometric pores and cavities. Indeed, it is well established that the properties of confined fluids is connected to a subtle interplay between the fluid-membrane interactions and geometrical restrictions.^1,2 These effects may drive a change in microscopic organization^3–5 and a drastic slowdown in translational and rotational dynamics.⁶ This last decade, a large body of theoretical and numerical research works has been devoted to investigate fluids under micro and nanoconfinement.

Molecular dynamics (MD) simulations is a numerical tool capable to describe water and ions at a molecular level and their behavior under confinement. There has been a large amount of works focused on electrolytic solutions trapped within nanotubes. While many works have provided molecular-level insights on the behavior of interfacial water,^7–11 on the physical properties of electrolyte solutions in bulk phase^12,13 and on aqueous solutions confined in hydrophobic frameworks^14–19 few works have been undertaken to study the structure and the dynamics of confined electrolyte solutions into hydrophilic frameworks.^4,20–23 Our motivation in performing this set of simulations was to explore the concentration dependence of both structural and dynamical properties of confined NaCl solutions at nanoscale into a hydrophilic pore model with a pore diameter of 1.2 nm. This simple curved geometry-model combines a radial restriction and an axial isotropic direction and allowed us to unravel the confinement and surface effects on the structure and dynamics of NaCl solutions.

The organization of the paper is a follow: computational details and force fields are presented in Sec. 2. The structural and dynamical properties of confined NaCl solutions are discussed in Sec. 3. Our concluding remarks are given in Sec. 4.

2 Computational details

Cylindrical silica nanotube was obtained by carving a hole of diameter D = 1.2 nm (see Fig. 1) in amorphous silica by means of the approach of Brodka and Zerda²⁴ leading to a hydrophilic nanotube with a density of surface silanol (SiOH) groups of 7.5 SiOH per nm². This matrix corresponds to a highly hydrated protonated silica pore.^24,25 Details on the framework building method and force field have been taken from elsewhere.^3,4,6,26,27 The axial length is L_z = 2L_x = 2L_y = 71 Å. The silica force field is provided in ref. 4. Although the silica matrix was subsequently kept rigid, rotation around the Si–O bond of the hydroxyl groups was allowed from the application of a bending potential between Si–O–H. Motion of hydrogen atoms of SiOH are performed from a stretching potential. Intramolecular parameters are given in ref. 28. Ions were modeled by means of an unpolarizable model²⁹ while the rigid TIP4P/2005 water model was considered.³⁰ Previously we showed that the polarizability of NaCl slightly impacts the structure and the dynamics of ions.^4,28 The density of confined water has been managed by means of a hybrid Grand Canonical Monte Carlo (GCMC)/Molecular dynamics (MD) simulations (GCMC/MD). Indeed, the concentration of electrolytes is kept constant. Therefore, to correctly model the dynamics of the ions and their hydration shell we implemented an MD move into our GCMC scheme. The theoretical background behind this hybrid approach is accurately detailed in ref. 31. For the MC procedure, the ratios for each trial move were defined as follows: 0.1995 for the translation, 0.1995 for the rotation, 0.6 for the insertion/deletion, and 0.001 for the MD move (corresponding to 2/N, usually considered as the frequency of the volume change). The MD runs corresponded to 5000 steps (N_MD) using a time step of 0.001 ps (t_MD). We considered more than 70 000 GCMC/MD cycles corresponding to 100 ns of simulation time. To validate our findings, several ratios for each trial move and (N_MD, t_MD) couples have been tested in order to confirm that the simulations do not depend on the initial conditions. The chemical potential (μ) of water molecules is 2431 K and corresponds to a pressure of 1 bar. Chemical potential was calculated from the Widom insertion during the simulation of the water bulk phase in the NpT statistical ensemble at 300 K and 1 bar.³² This method consists in inserting a ghost particle randomly into a simulation box of N molecules and calculating the interaction energy of this ghost particle with the N molecules. Molecular dynamics (MD) simulations were performed using a time step of 0.001 ps to sample 10 ns (acquisition phase). All MD simulations have been carried out with the DL_POLY package³³ using a combination of the velocity-Verlet algorithm and the Nose–Hoover thermostat^34,35 with a relaxation time of τ_t = 0.5 ps. Periodic boundary conditions were applied in the three directions. MD simulations were performed in the NVT statistical ensemble at T = 300 K. Molecular dynamics (MD) simulations were performed using a time step of 0.001 ps to sample 10 ns (acquisition phase). The equilibration time corresponds to 10 ns. The equation of motions was integrated using the velocity Verlet scheme coupled with the QUATERNION, SHAKE-RATTLE algorithms.³⁶ The Ewald summation was used for calculating the electrostatic interactions and the short range interactions were truncated at 12 Å. The statistical errors for the calculated properties (as the diffusion coefficient) are estimated using the block average method. This method uses four superblocks which are formed by combining three blocks of 2000 configurations. The radial profiles were calculated by using 4 blocks of configurations. The error bars on radial profiles are too smalls to be indicated. The comparison between the bulk and the confined properties was performed at T = 300 K and p = 1 bar. Six concentrations in NaCl were investigated {0, 0.2, 0.4, 1.0, 2.0, 5.0} mol L⁻¹.

Fig. 1 Axial view of the final configuration of the confined NaCl solution at 5 mol L⁻¹. Red and yellow colors correspond to the oxygen and silicon atoms, respectively. Hydrogen atoms are in white. sodium and chloride ions are in blue and green respectively.

3 Results and discussion

3.1 Structural properties

The structure of confined electrolyte solutions was first investigated by considering the distribution of the species inside the pore. Fig. 2 shows the radial profile of water density inside the pore for different salt concentrations. Two maxima are evidenced, which suggests a layering organization into the pore induced by the confinement effect.¹¹ As shown in Fig. 2 three boundaries can be drawn: two zones corresponding to both shells (the center of pore between 0 and 2 Å and the interfacial region between 3 and 6 Å) and a third region between the center of pore and the interfacial region between 2 and 3 Å corresponding to the interlayer domain. From c = 0 mol L⁻¹ (pure water) to c = 2 mol L⁻¹ the first peak corresponding to the center of pore is shifted toward the silica surface while the second peak located at 4.8 Å is unaffected by the increase in NaCl concentration. On the other hand, a shift of both layers towards the pore surface in observed for c = 5 mol L⁻¹. Interestingly, let us note the intensity of the second peak remains unchanged despite a change in water molecules number into the pore as function of the NaCl concentration (obtained from GCMC simulations). This suggests that the interfacial density is ruled by the surface–water interactions whereas water density in the center of pore decreases as NaCl concentration increases. Additionally, we clearly show in Fig. S1a and S1b of the ESI† that both Na⁺ and Cl⁻ ions are located at the center of the pore, i.e. in the first water layer. This explains why water density in the pore center is strongly impacted by the addition of salt whereas the pore surface is almost not altered. Furthermore, the shift of the location of the maximum density of water for c = 5 mol L⁻¹ can be correlated to the local maximum of the ion density observed at 3 Å. In Fig. S1† we show a slight difference between the local densities of Na⁺ and Cl⁻ ions. Despite this small difference in the local density, we checked through the radial distribution in Fig. S4a† the ions pairs formation. Let us note that the local concentration was calculated by calculating the density distribution into a cylindrical shell of radius of 0.2 Å.

Fig. 2 Radial water density profile. The pore center is located at r = 0 Å and the pore surface is around r = 6 Å. We defined three regions, the centre of pore, the interlayer domain and the interfacial region.

To investigate the local structure of water molecules the distribution of the tetrahedral order parameter () was computed (P(q)).³⁷ ϕ_ij is the angle between i and j molecules. For a tetrahedral structure q is close to 1. In Fig. 3 we report P(q) at the pore center and close to the pore surface. As shown in Fig. 3a the local organization of water at the pore center is similar to that in the bulk phase (Fig. S2 of the ESI†), whatever the salt concentration. At the pore center and into the bulk phase the presence of ions modify significantly the local tetrahedral organization of water molecules as shown by the shift of the most probable value of q. Indeed, the most probable value of q is shifted from q = 0.8 for pure water to around 0.5 for c = 5 mol L⁻¹. The modification of the tetrahedral organization of water molecules is imputed to the strong orientation of water molecules around ions in both bulk phase and the central part of the pore. As shown Fig. 3b a different behavior is observed close to the pore surface and the most probable value of q is found to be 0.5, whatever the concentration. This local structure is due to water–surface interactions which orientate water molecules (see discussion below). To sum up, the local organization of water molecules is mainly impacted by the pore surface in the interfacial region while it is disturbed by ions in the central part of the pore.

Fig. 3 Tetrahedral order distribution of water molecules at the pore center (a) and in the interfacial region (b) as function of the NaCl concentration. These two regions are defined in Fig. 2.

In order to investigate the local structure around ions and water molecules radial distribution functions (RDF) were computed. All RDF were corrected with the excluded volume.¹ Fig. S3a and S3b† show the radial distribution functions of water oxygen atoms around Na⁺ ions and water hydrogen atoms around Cl⁻ ions, respectively. Whatever the NaCl concentration location of the maxima is the same, which suggests that the nature of the interactions is kept. In Fig. S3c and S3d of the ESI† the profile of hydration number of ions as function of the radial position is reported. These numbers were calculated from the integration of the radial distributions functions presented in Fig. S3a and S3b.† In the central part of the pore, the coordination number reaches the bulk value for both ions i.e. 5.7 and 6.5 in the first hydration shell for Na⁺ and Cl⁻, respectively, for concentrations up to 2 mol L⁻¹.² For the highest concentration (c = 5 mol L⁻¹) a strong decrease is observed at the pore center with respect to the bulk value. In this case none of the ions is bulk-like hydrated. This is due to the confinement effect which favors the formation of Na–Cl clusters. As shown in Fig. S3c and S3d of the ESI† the hydration number of ions decreases approaching the surface. This is due to the presence of the surface that truncates the hydration shell of ions. At 5 mol L⁻¹ the hydration number of Na⁺ is constant along the radial position (r). This is the consequence of the formation of the ions aggregates at the centre of pore which decrease the hydration number from 6 to 4 and which corresponds to the hydration number of ions close to the surface. To investigate ion-pairing we report in Fig. S4a of the ESI† the RDF between Na⁺ and Cl⁻. The first shell is located at r = 3.5 Å whatever the concentration. The same distance is obtained in bulk phase, which suggests a weak impact of confinement on the intrinsic nature of interactions between ions. RDFs between Na⁺ ions and their integration are reported in Fig. S4b of the ESI† which highlights the formation of ion clusters into the pore. The maximum of the RDF around r = 4 Å corresponds to the situation where one chloride anion is surrounded by 2 sodium cations. In addition, the local maxima observed in radial density in ions in Fig. S1 of the ESI† can be imputed to a formation of cluster ions.

The analysis of hydrogen bonds (HB) between two water molecules and between a water molecule and the silica surface was performed by means of a geometrical criterion based on quantum calculations.³⁸ Two molecules are bonded if (i) the oxygen–oxygen distance is less than 3.5 Å, (ii) the distance H–O is less than 2.5 Å and (iii) the oxygen–oxygen–hydrogen angle is less than 30°. The exploration of interactions between water molecules and ions was performed by using another criterion.¹³ A chloride anion and a water molecule are considered bonded if (i) the oxygen–chloride distance is less than 3.90 Å, (ii) the hydrogen–chloride distance is less than 3.05 Å and (iii) the chloride–oxygen–hydrogen angle is less than 45°. Fig. 4 shows the radial profile of the total number of hydrogen bonds per water molecule (nHB). This calculation includes 3 types of HBs water–water, water–silica and water–ions. In general nHB decreases as the concentration increases. At the center of pore (between 0 and 2 Å) this decrease is imputed to the ions which orientate water molecules around them and then disturb the HB network by breaking water–water HBs. This decrease in water–water HBs is not compensated by the physical bonds between ions and water molecules (chaotropic effect). As shown in Fig. 4 the decrease in nHB as function of the NaCl concentration is aggravated under confinement with respect to the bulk phase due to the excluded volume effect which prevents to recover the bulk HB coordination. In the interlayer zone an unexpected increase in nHB is observed. This rise can be imputed to the additional HBs between the two water layers. In the interfacial zone the number of HB decreases as function of r due the excluded volume effect and reaches the same value, whatever the salt concentration. The interfacial value is then set by the SiOH surface. At 5 mol L⁻¹ the HB profile is constant in the interfacial region because the limiting interfacial nHB value is already reaches from r = 3.5 Å. This structural behavior highlights the subtil interplay between surface and confinement effects and the presence of ions inside the confining medium. The radial profile of the HB energy is shown in Fig. S5 of the ESI.† In this calculation the HB between two water molecules was considered. We observe that the HB energy of confined water is lower than the one of the bulk phase. The energetic difference between the bulk and the confined phases increases with the NaCl concentration, which suggests a decrease in the cohesion of the HB network under confinement.

Fig. 4 Radial profile of the total number of hydrogen bonds per water molecule as function of the NaCl concentration. Bulk values are shown by the horizontal lines.

3.2 Dynamical properties

The translational dynamics was investigated by computing the self-diffusion coefficients of each species inside the pore. The Einstein's relation was used to compute the self-diffusion coefficient from the mean square displacement (MSD). The reader is directed to ref. 39–41 for some details on the self diffusion coefficient calculation. As underlined by Lindahl and Edholm⁴⁰ and by Sega et al.⁴¹ calculation of diffusion coefficient in confined regions is difficult given the non-linearity of the MSD at longtime (MSD ∼ Dt^α with α ≠ 1).³⁹ During the last 2 ns we found α = 1.018 for the confined water. This result allows us to use the Einstein's relation for a linear regimewith t₀ the time origin, N the number of molecules, N₀ the number of t₀ and r_com,i the position vector of the center of mass of i molecule. n is the dimension of the system such as n = 3 in the 3D isotropic phase. Under cylindrical confinement n = 1 and n = 2 in the axial and radial directions, respectively. Fig. 5 shows the MSD of water molecules as a function of NaCl concentration inside the pore. First of all we observe a dramatic slowdown of water translational dynamics between bulk (Fig. S6 of the ESI†) and confined phases. In the bulk phase the MSD of water varies linearly with time whatever the NaCl concentration while inside the pore a crossover is observed (at t = 9 ps for pure water and up to t = 40 ps for the 5 mol L⁻¹ solution). This latter corresponds to a change in the translational dynamics from a subdiffusive regime at short times to a fully diffusive regime at longer times (MSD ∼ t) which is a typical feature of dense liquids. Indeed, as shown in Fig. 2 the confined water density increases in relation to the bulk phase. The subdiffusive regime is more pronounced as the NaCl concentration increases. Indeed, the subdiffusive regime is associated to the cage effect where a water molecule is locked up in a water cage while the diffusive regime corresponds to the release of water molecule. In the electrolytic solutions Na⁺ and Cl⁻ ions hinder the diffusion of water molecules and delay the escape from the cage.

Fig. 5 Mean square displacement (MSD) of water molecules in confined phase as function of the NaCl concentration. For 5 mol L⁻¹, the two solid lines are guide lines to highlight the crossover between the subdiffusive and the diffusive regime.

The self-diffusion coefficients were extracted from the last 2 ns of the MSD i.e. during the diffusive regime. We report in Fig. 6 the self-diffusion coefficient of water in bulk and confined phases as function of the salt concentration. In both cases we show that the presence of ions slows down the translational diffusion of water molecules, this slowdown being all the more pronounced as the NaCl concentration increases. This effect is due to the strong orientation of water molecules around ions leading a decrease in the translational dynamics of water molecules. Let us note that the self-diffusion coefficients of ions were not extracted because the statistical fluctuations were too high given the small number of ions inside the pore.

Fig. 6 Water self-diffusion coefficients as function of the NaCl concentration in the bulk and confined phases. Diffusion coefficient was calculated by means of eqn (1) with n = 3 and n = 1 for the bulk and confined phases, respectively.

The rotational dynamics of the water molecules was investigated from the calculation of the dipolar relaxation time (τ_d) of water obtained by means of an adjustment of the correlation function of the total dipole moment.⁴ Two areas in the confined phase were distinguished, one at the pore center and another one close to the surface. Relaxation times are reported in Table 1.

Table 1 Dipolar relaxation time (τ) (ps) and diffusion coefficient (D) (10⁻⁹ m² s⁻¹) of water in the bulk, at the center of pore and close to the surface. Results are given in ps. Error bars of the diffusion coefficients and the relaxation times are about 0.1 × 10⁻⁹ m² s⁻¹ and 0.2 ps respectively

	Bulk	Pore center	Pore surface
	τ
Pure water	4.69	6.37	8.86
0.2 mol L^−1	4.85	8.58	11.57
0.4 mol L^−1	5.09	9.23	12.40
1 mol L^−1	5.63	11.24	13.46
2 mol L^−1	7.05	15.49	16.35
5 mol L^−1	12.8	89.27	99.99

---

	Bulk	Pore center	Pore surface
	D
Pure water	2.39	0.31	0.31
0.2 mol L^−1	2.30	0.30	0.30
0.4 mol L^−1	2.21	0.29	0.29
1 mol L^−1	1.95	0.26	0.26
2 mol L^−1	1.54	0.22	0.25
5 mol L^−1	0.81	0.09	0.1

---

Due to confinement effect, τ_d of confined water molecules is higher than in the bulk phase. Moreover, the relaxation time increases as the ions concentration due to the local dipolar orientation of water in the vicinity of ions. As in the bulk phase, at the center of pore the ions slow down the rotational motion of water and this effect is stressed for the high concentrations. Indeed, at c = 5 mol L⁻¹ the relaxation time of confined water is 89.3 ps while in the bulk phase τ_d = 12.8 ps. This increase in τ_d is due to a combination of the salt effect and confined environment which induces a slow down in the rotational motion. Close to the pore surface the relaxation time of water molecules is of the same order of magnitude as at the center of pore despite the absence of ions in the vicinity of the surface. This result can be explained by the SiOH surface which induces a strong hydrophilic anchoring that slows down the rotational motion of water molecules. Dipolar orientation of water molecules can be studied by means of the distribution of the angle (Θ) formed by the dipole moment vector and the vector normal to the silica surface. We report in Fig. 7 the distribution of Θ at the pore center (between 0 and 1 Å) and close to the surface. At the pore center a flat profile is observed while near to the interfacial a maximum is located around 98°. This preferential orientation results from water–SiOH interactions which orientate water molecules to favor HB formation. It can be mentioned that this preferential orientation propagates from the surface to r = 2 Å (not shown). This is borne out by the similar dipolar relaxation times and diffusion coefficients computed at the pore center and close to the surface (Table 1). Interestingly, we note in Table 1 the strong gap between the dipolar relaxation time of water at c = 2 mol L⁻¹ and c = 5 mol L⁻¹. This difference is directly imputed to Na–Cl clusters which strongly orientate water molecules.

Fig. 7 Distribution of Θ as function of the NaCl concentration (a) at the center of pore and (b) in the interfacial region.

4 Concluding remarks

Structure and dynamics of NaCl solutions confined in a hydrophilic cylindrical (diameter: 1.2 nm) nanotube were investigated by means of molecular dynamics simulations. We showed that the modifications of the structural and dynamical properties of confined electrolyte solutions result from the interplay of confinement, surface and ion effects. The geometric restrictions (confinement effect) lead to the structuration of water molecules in two concentric layers. Water density in the interfacial layer was found to be imposed by surface–water interactions while water density in the central part of the pore was mainly affected by the presence of ions.

The presence of ions modifies significantly the local tetrahedral organization of water molecules. Close to the pore surface, the local organization of water molecules is mainly impacted by the pore surface while ions affect the local order of water molecules in the central part of the pore. The decrease in hydrogen bonds number as function of the NaCl concentration is enhanced under confinement with respect to the bulk phase. This structural behaviors highlights the subtle interplay of surface, confinement and salt effects. We also showed that ion pairs form at the pore center so as to keep their hydration shells and that confinement favors the formation of Na–Cl clusters with respect to the bulk phase.

From an energetic standpoint, the energy of water–water HB is weaker inside the nanotube than in the bulk phase. This energetic difference between bulk and confined phases increases with NaCl concentration, which suggests a decrease in the cohesion of the HB network under confinement.

Confinement decreases significantly water diffusivity and a subdiffusive regime was evidenced inside the nanotube. This subdiffusive regime occurs over longer times as the NaCl concentration increases because Na⁺ and Cl⁻ ions hinder the diffusion of water molecules and delay the escape from the cage. The presence of ions slows down the translational and rotational dynamics of water molecules in both bulk and confined phases. Finally, confinement was found to strengthen the impact of ions on rotational dynamics of water molecules.

Acknowledgements

The authors are grateful to the Agence Nationale de la Recherche for its financial support through the program MUTINA (ANR 2011 BS09 002).

References

1. D. Morineau and C. Alba-Simionesco, J. Chem. Phys., 2003, 118, 9389–9400.
2. C. A. Simionesco, B. Coasne, G. Dosseh, G. Dudziak, K. Gubbins, R. Radhakrishan and M. Sliwinska-Bartowiak, J. Phys.: Condens. Matter, 2006, 18, R15.
3. A. Ghoufi, I. Hureau, R. Lefort and D. Morineau, J. Phys. Chem. C, 2011, 115, 17761.
4. R. Renou, A. Ghoufi, A. Szymczyk, H. Zhu, J.-C. Neyt and P. Malfreyt, J. Phys. Chem. C, 2013, 117, 11017.
5. A. Ghoufi, I. Hureau, D. Morineau, R. Renou and A. Szymczyk, J. Phys. Chem. C, 2013, 117, 15203.
6. A. Ghoufi, D. Morineau and R. Lefort, J. Chem. Theory Comput., 2010, 6, 3212.
7. P. Gallo, M. Rapinesi and M. Rovere, J. Chem. Phys., 2002, 117, 369–375.
8. N. Giovambattista, P. G. Debenedetti and P. J. Rossky, J. Phys. Chem. B, 2007, 111, 9581–9587.
9. J. Jelassi, T. Grosz, I. Bako, M.-C. Bellissent-Funel, J. C. Dore, H. L. Castricum and R. Sridi-Dorbez, J. Chem. Phys., 2011, 134, 064509.
10. A. A. Milischuk and B. M. Ladanyi, J. Chem. Phys., 2011, 135, 174709.
11. R. Renou, A. Szymczyk and A. Ghoufi, J. Chem. Phys., 2014, 140, 044704.
12. R. D. Mountain, Int. J. Thermophys., 2007, 28, 536–543.
13. A. Nag, D. Chakraborty and A. Chandra, J. Chem. Sci., 2008, 120, 71–77.
14. B. Corry, J. Phys. Chem. B, 2008, 112, 1427.
15. T.-A. Beu, J. Chem. Phys., 2011, 135, 044516.
16. M. Wander and K. Shuford, J. Phys. Chem. C, 2010, 114, 20539.
17. P.-A. Cazade, J. Dweik, B. Coasne, F. Henn and J. Palmeri, J. Phys. Chem. C, 2010, 114, 12245–12257.
18. L. Wang, R. S. Dumont and J. M. Dickson, J. Chem. Phys., 2012, 137, 044102.
19. L. Wang, R. S. Dumont and J. M. Dickson, J. Chem. Phys., 2013, 138, 124701.
20. P. E. Videla, J. Sala, J. Martí, E. Guàrdia and D. Laria, J. Chem. Phys., 2011, 135, 104503.
21. J. Dweik, B. Coasne, J. Palmeri, P. Jouanna and P. Gouze, J. Phys. Chem. C, 2012, 116, 726–737.
22. H. Kulik, E. Schwegler and G. Galli, J. Phys. Chem. Lett., 2012, 3, 2653.
23. P. A. Bonnaud, B. Coasne and R. J.-M. Pellenq, J. Chem. Phys., 2012, 137, 064706.
24. A. Brodka and T. Zerda, J. Chem. Phys., 1996, 104, 6319.
25. J. Puibasset and J.-M. R. Pellenq, J. Chem. Phys., 2005, 122, 094704.
26. R. Guégan, D. Morineau and C. A. Simionesco, J. Chem. Phys., 2005, 317, 236.
27. R. Busselez, R. Lefort, Q. Ji, F. Affouard and D. Morineau, Phys. Chem. Chem. Phys., 2009, 11, 11127.
28. H. Zhu, A. Ghoufi, A. Szymczyk, B. Belannec and D. Morineau, Phys. Rev. Lett., 2012, 109, 107801.
29. R. G. Moutain, Int. J. Thermophys., 2007, 28, 536.
30. J. Abascal and C. Vega, J. Chem. Phys., 2005, 123, 23505.
31. A. Ghoufi and G. Maurin, J. Phys. Chem. C, 2010, 114, 6496.
32. A. Ghoufi, D. Morineau, R. Lefort, I. Hureau, L. Hennous, H. Zhu, A. Szymczyk, P. Malfreyt and G. Maurin, J. Chem. Phys., 2011, 134, 074104.
33. I. Todorov, W. Smith, K. Trachenko and M. Dove, J. Mater. Chem., 2006, 16, 1911.
34. S. Nose, J. Chem. Phys., 1984, 81, 511.
35. W. Hoover, Phys. Rev. A, 1985, 31, 1695.
36. M. P. Allen and D. Tildesley, Computer simulations of liquids, Oxford, 1987.
37. J. Errington and P. Debenedetti, Nature, 2001, 409, 318.
38. A. Luzar and D. Chandler, J. Chem. Phys., 1993, 98, 8160–8173.
39. A. Bizzarri and S. Cannistraro, Phys. Rev. E, 1996, 53, R3040.
40. E. Lindahl and O. Edholm, Phys. Rev. E, 1998, 57, 791.
41. M. Sega, R. Vallauri and S. Melchionna, Phys. Rev. E, 2005, 72, 041201.

---

Footnote

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

---