Complex Behavior of Ordered and Icelike Water in Carbon Nanotubes

Letter pubs.acs.org/JPCL

Cite This: J. Phys. Chem. Lett. 2018, 9, 4746−4752

Complex Behavior of Ordered and Icelike Water in Carbon Nanotubes near Its Bulk Boiling Point Jose ́ Cobeña-Reyes, Rajiv K. Kalia, and Muhammad Sahimi* Mork Family Department of Chemical Engineering and Materials Science, University of Southern California, Los Angeles, California 90089-1211, United States

Downloaded via KAOHSIUNG MEDICAL UNIV on September 3, 2018 at 06:09:36 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: We report the results of extensive molecular dynamics (MD) simulation of water in a carbon nanotube (CNT) with a specific diameter over a wide range of temperatures from 343 to 423 K. In order to characterize the nature of water, we have computed the Kirkwood g-factor, the ten Wolde parameter, the radial distribution, the cage correlation, the intermediate scattering functions, the mean-square displacements of the water molecules, and the connectivity of the oxygen atoms. The computed properties provide evidence for complex behavior. Some of the properties indicate an icelike structure, while others point to ordered (but not necessarily frozen) water. The connectivity is close to 9. The ordered water exists both below and above its bulk boiling point. The order is identified based on the ten Wolde parameter and may explain, along with the dynamic slow down, the recent discovery of “ice” in CNTs near the bulk boiling point in a certain range of CNT diameters, not seen in tubes of other sizes.

t is well-known that some fluids exhibit anomalous behavior, which is different from both experimental observations and theoretical expectations for most liquids. For example, whereas the diffusion coefficient of most liquids decreases under compression, there are some fluids whose diffusivity increases under the same conditions. Examples include, but are not limited to, Si,2,3 graphite,4 and liquid metals.5 The most important fluid with anomalous behavior is, however, water,6,7 whose density is maximum at 277 K; it may possibly coexist within two distinct forms8 (see, however, Selberg et al.9 and Laksmono et al.10); it is in a supercooled state at sufficiently low temperatures,8 and the behavior of its isothermal compressibility, isobaric heat capacity, and thermal expansion coefficient11 at temperatures between homogeneous nucleation at 231 K and its melting point at 273 K is unlike that of almost any other liquid. Moreover, it is well-known that three forms of glassy water exist under bulk conditions, namely, low-density, high-density, and very high density amorphous ice, and that, unlike most other liquids, volumetric expansion of water at 277 K and lower temperatures is due to the reduction in its entropy caused by the tetrahedral symmetry of the local order around each water molecule12 and the formation of hydrogen bonds. Given that such anomalous behavior of water occurs in an unbounded (bulk) medium, an important question to be addressed is what are water’s properties in confined media and, in particular, in nanotubes. One cannot overstate the importance of the question because understanding the properties of water in nanostructured materials and media is not only important from a fundamental theoretical viewpoint but is also critical to many physical, chemical, and biological phenomena. Experiments11,12 as well as computer simulations indicate that the water flux in carbon nanotubes13 (CNTs) and

I

© 2018 American Chemical Society

silicon carbide nanotubes14 (SiCNTs) may be 3−4 orders of magnitude larger than what is predicted by classical Hagen− Poiseuille flow (HP) in cylindrical tubes. In addition, it has been demonstrated that other factors, such as functionalization of nanotubes, strongly influence water flow in CNTs.15,16 Similarly, whereas diffusion of water in unbounded media follows the Einstein relation and is affected by only the local temperature and pressure, in a confined medium, the van der Waals and Coulombic interactions between water molecules and the medium’s walls influence their mobility, which typically reduces the rates of diffusion.17,18 The Stokes− Einstein relation between the diffusivity and viscosity of water also breaks down in confined media,19 including CNTs and SiCNTs.20,21 Using Raman spectroscopy, Agrawal et al.22 studied the behavior of water in isolated CNTs. They reported very high sensitivity to the diameter of CNTs and much larger increases in the freezing transition temperature than what has been predicted theoretically. Using six isolated CNTs of various diameters between 1.05 and 1.52 nm, Agrawal et al.22 reported various unusual phenomena that depended on the size of CNTs, including reversible melting bracketed to 105−151 and 87−117 °C, near-ambient phase transitions bracketed between 15 and 49 and 3−30 °C, and depression of the freezing point between −35 and 10 °C. Such unusual behavior had never been reported before. In this Letter, we report the results of extensive molecular dynamics (MD) simulation of water in CNTs over a wide Received: June 21, 2018 Accepted: August 3, 2018 Published: August 3, 2018 4746

DOI: 10.1021/acs.jpclett.8b01953 J. Phys. Chem. Lett. 2018, 9, 4746−4752

Letter

The Journal of Physical Chemistry Letters range of temperature and with the same size that Agrawal et al.22 studied. We present the results for several important properties of water in CNTs, which provide evidence for a highly complex behavior of water, including the presence of ordered water with connectivity much higher than its connectivity under bulk conditions and icelike structures that are similar to the purported ice reported by Agrawal et al.22 We have computed several quantities in order to characterize the state of water in the CNTs. In what follows, we present and discuss the results. Kirkwood g-Factor. We computed the Kirkwood g-factor for two purposes. One was, as described in the Supporting Information (SI), to determine the appropriate size L of the simulation cell and, hence, that of CNTs l in order to produce results that are independent of l. The results for this purpose are presented in the SI. The second purpose was to gain insight into the local orientational order of water that the g-factor measures. Thus, we computed the change of the g-factor as the temperature was lowered from 423 K. Figure 1 presents the

Figure 2. Comparison of the ten Wolde parameter for the bulk and confined water at T = 363 K.

neighbors present in liquid water under bulk conditions. Thus, one has a mixture of liquid and icelike states. Dependence of the Potential Energy on the Radial Distance. An interesting feature is the dependence of the potential energy U(r) on the radial distance r. As Figure S3 in the SI indicates, there are essentially two groups of water molecules in CNTs. One group is roughly in a cylinder near the CNT’s center, while the second group is near the wall. Thus, we split the water molecules into concentric cylinders every 0.5 Å in the radial direction, starting at the center, r = 0, and computed the potential energy U(r) of each group. As Figure 3 indicates, the

Figure 1. Kirkwood g-factor at several temperatures.

results. The g-factor at high temperatures T indicates no particular dipole orientation (see also Figure S2 in the SI). As T decreases, however, the g-factor begins to increase, indicating that there is a dipole orientation parallel to the van der Waals force and the interaction between the CNT and the water molecules. Because the water molecules are in a lowerfrequency mode (they vibrate less) at lower temperatures, such interactions become relatively stronger, producing an increase in the g-factor and the icelike structure. Note the gap between the g-factor for the three lowest temperatures and the higher temperatures that begin at at 373 K, indicating different water structure at the lower temperatures. Note also that the highest g-factor is at 363 K. The increase in the g-factor is consistent with the increase in the ten Wolde parameter and the slow down in the mean-square displacements (MSDs), to be described shortly. ten Wolde Parameter. As explained in the Methods section, we computed the ten Wolde parameter dl(i,j) for l = 6. Figure 2 compares the results for water in CNTs at two temperatures with the corresponding results for bulk water. There are very significant differences between the confined water order parameter and that in the bulk. Note that values of d6(i,j) for about 50% of the water molecules are larger than 0.50, implying that at least about 50% of the molecules are coherently “bonded.” In addition, at least 70% of the water molecules have 8−9 neighbors, much larger than the typical 4

Figure 3. Dependence of the potential energy U(r) on the radial distance r from the CNT’s center.

water molecules near the CNT’s wall possess much lower values of U(r) than those near the center, indicating a slow down in the dynamics in that layer. Such a difference in the dynamics was also confirmed by calculating the cage correlation (CC) function, which will be discussed next. Cage Correlation Function. Figure 4a presents the computed CC functions C(t) at four temperatures. Although C(t) decays with time, even after 40 ns it is still significantly larger than zero. Similar to the potential energy, we split the CNT into two cylinders. Any water molecule between r = 0 (the center) and 1.25 Å was considered to belong to the cylinder around the tube’s center. All other water molecules belonged to what we refer to as the outer layer, which includes those near the wall. Figure 4b presents the results for the two groups. The CC function for the outer layer decays much slower than that near the center, indicating that any motion to break the cage in the outer layer is far less likely, or much slower, than those near the 4747

DOI: 10.1021/acs.jpclett.8b01953 J. Phys. Chem. Lett. 2018, 9, 4746−4752

Letter

The Journal of Physical Chemistry Letters

Figure 4. CC function C(t). (a) Overall function at several temperatures. (b) Comparison of C(t) for the water molecules near the tube’s center and in the outer layer.

exponential decay, characterized by a relaxation time τl. This is a typical glasslike behavior where the Gaussian feature vanishes, making the stretched exponential decay dominant as time increases, which happens for large κ. In this case, at higher temperatures, the function decays faster, but again, even after 40 ns, it has not vanished, except for the values that correspond to short distances. Fitting the data to eq 14 yielded τl ≈ 88 ns for κ = 2, the highest peak in the structure factor. Mean-Square Displacements. The computed results for the axial MSDs indicate a drastic slow down in the dynamics of the water. Indeed, if we consider the MSDs shown in Figure 6, we

center. In other words, the dynamics only somewhat away from the center is much slower than that near it. Because the ten Wolde parameter indicates the presence of ordered structures, with 50% or more of the water molecules exhibiting icelike behavior and the molecules in the outer layer also having very low potential energies than those near the center, as well as a very slowly decaying CC function, we may conclude that an icelike phase exists in that region. Thus, two phases of water, icelike and liquid, may coexist in the CNT, consistent with the observation that water coexists in different phases in a CNT.1 Intermediate Scattering Function. Figure 5 presents Fs(κ,t) as a function of time t and κ. As expected, the decay is faster with

Figure 6. Axial MSDs ⟨R2(t)⟩ of water (oxygen) molecules. Figure 5. Intermediate scattering function at 343 K.

see that at T = 343 K they grow with time very slowly and after some time saturate. A similar phenomenon happens at T = 353 K, albeit after a somewhat longer time. The maximum values of the MSDs at the two temperature are around 10−12 Å2, implying a net displacement of about 3.3−3.5 Å, which is about the only Lennard-Jones (LJ) size parameter used in various models of water and essentially represents vibration of the water molecules around their positions and not truly significant displacements. However, at higher temperatures, the MSDs appear to rise with time without any saturations over the time period in which we carried out our simulations. Radial Distribution Function and Mean Connectivity of the Atoms. Figure 7 presents the computed radial distribution functions (RDF) at T = 363 K for the pairs O−O, H−H, and H−O. The results for T = 343 and 353 K are completely similar. While the general trends in all three RDFs are similar, there are also subtle differences. The main peak in g(r)OO at

high values of κ (short distances) but is slow at low wavenumbers κ (long distances). Qualitatively, the data shown in Figure 5 are similar to those reported in several studies of supercooled or glasslike water23−26 and confined water at very low temperatures27,28 in various materials. The shape of Fs(κ,t) and the relaxation times that we have computed are also comparable to those in the aforementioned studies. As Figure 5 indicates, Fs(κ,t) behaves qualitatively similar to the CC function C(t) in that it indicates a notable slowing down of the dynamics. That is, the water molecules go thorough two stages of relaxation. The first one is a fast process that is described by a Gaussian decay in which the molecules move inside of the cage produced by the other molecules,25 whereas the second relaxation is described by a slow stretched 4748

DOI: 10.1021/acs.jpclett.8b01953 J. Phys. Chem. Lett. 2018, 9, 4746−4752

Letter

The Journal of Physical Chemistry Letters

experiments of Agrawal et al.22 The CNT was then inserted inside of a simulation cell of size 3502 × l Å3, filled with enough water molecules to have a water density of 1 g/cm3. We then proceeded with energy minimization, followed by thermalization. The minimization was done in three steps, using the conjugate-gradient method. First, we minimized only the energy of the water molecules, followed by that of the CNT only, and finally for the entire system in order to prevent the atoms from moving large distances over short times. In the thermalization step, we increased the temperature using the Langevin thermostat and carried out the simulations in the (NVE) ensemble to reach the equilibrium state. The temperature increase for all of the systems that we simulated was 1 K after every 2 ps. The total time for energy minimization and thermalization was about 1.2 ns. After reaching the desired temperature, MD simulations were carried out using the velocity-Verlet algorithm in an (NVT) ensemble for 40 ns. All of the simulations were carried out using the LAMMPS package,29 and the system was visualized using Chimera.30 The CNT contained 3620 carbon atoms, and there were 851 water molecules in the nanotube. The time step was always 1 fs. The REBO potential31 was used to model the CNT, while the TIP4P-EW32 was utilized to describe the water molecules. The latter was used because the water boiling point that it predicts, 370.3 ± 1.9 K, is the closest to the experimental bulk value.32 The van der Waals interactions between the carbon atoms and the water molecules were represented by the LJ potential with a cutoff distance of 11 Å. The LJ parameters for graphite,33 σ = 3.4 Å and ϵ/kB = 28.1 K−1, were utilized for the carbon atoms, where kB is the Boltzmann factor. The Lorentz− Berthelot mixing rules, σij = (σii + σjj)/2 and ϵij = ϵiiϵjj , were utilized for the pairs ij. The Coulombic interactions were computed using the particle−particle−particle−mesh method. Computing the Kirkwood g-Factor. When simulating any phenomenon in strong confinement, it is important to carry out the computations in such a way that any possible artifact is eliminated or at least minimized. Because we study water in nanotubes, one possible artifact is an artificial crystallinity that may arise due to the periodicity and summations methods for long-range interactions, such as the Ewald or particle− particle−particle−mesh method.34 Previous studies35−38 showed that the Kirkwood g-factor gK is an accurate tool for checking the existence of this type of artifact. gK is a measure of the orientational order of the dipole moment of the molecules and represents the correlation of the dipole μi of water molecule i with the total dipole of the water molecules located within a distance r from the molecule.37 gK = 1 when the dipole moment of the molecules lacks a specific orientation; for gK > 1, the dipole moments are parallel to an imposed field, whereas gK < 1 represents a case in which the dipole moments are aligned perpendicular to the field.39 In addition to gaining insight into the problem that we study, we also used gK to identify the right size of the simulation box. The g-factor is defined as39,40 μi ·μj gK (r ) = ∑ 2 μ r