THE JOURNAL OF CHEMICAL PHYSICS 142, 134704 (2015)
A transition between bistable ice when coupling electric field and nanoconfinement Feng Mei,1,2 Xiaoyan Zhou,1,3 Jianlong Kou,1 Fengmin Wu,1,3 Chunlei Wang,4,a) and Hangjun Lu1,a) 1
Department of Physics, Zhejiang Normal University, Jinhua 321004, China Department of Physics, Zhejiang University, Hangzhou 310027, China 3 Department of Physics and Institute of Theoretical Physics, Shanxi University, Taiyuan 030006, China 4 Division of Interfacial Water and Key Laboratory of Interfacial Physics and Technology, Shanghai Institute of Applied Physics, Chinese Academy of Sciences, P.O. Box 800-204, Shanghai 201800, China 2
(Received 16 January 2015; accepted 17 March 2015; published online 2 April 2015) The effects of an electric field on the phase behavior of water confined inside a nanoscale space were studied using molecular dynamics simulations. It was found that the diffusion coefficient of water reaches its maximum when value of the surfaces’ charge is at the threshold, qc = 0.5e. This unexpected phenomenon was attributed to the intermediate state between two stable ice states induced by nanoconfinement and the electric field generated by charged surfaces, respectively. Our finding is helpful to understand electromelting and electrofreezing of water under nanoconfinement with the electric field. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4916521]
I. INTRODUCTION
Understanding, controlling, and predicting the crystallization process of nano-confined liquids are fundamental to physics, chemistry, biology, and engineering applications.1–5 Controlling crystallization in cells can determine the ability of many organisms such as polar fish, insects, bacteria, and plants, to survive and thrive in subfreezing habitats.6 In the food, pharmaceutical and chemical industries, controlling crystallization in confined spaces is critical for both functional performance and regulatory compliance.4 However, the crystallization behavior of water is a complex phase transition that presents a challenge to statistical mechanics. Confinement of water inside nanoscale geometries or in the vicinity of surfaces can lead to interesting properties and increase the richness of its phase behavior, which differs considerably from bulk systems.7–26 A number of additional ice phases have been observed by simulating water confined within carbon nanotubes.7,13,16,18,27–32 Experimental results showed that decreasing the pore diameter depressed the freezing/melting point according to a modified Gibbs–Thomson equation.33,34 Computational and experimental studies have demonstrated the spontaneous freezing of monolayer, bilayer, and three-layer ice confined between parallel plates separated by nanometer gaps.7,26,35–38 Interestingly, there is still a debate about whether the nanoscale confinement can favor or suppress crystallization. Another factor that may affect water crystallization is surface polarity, which is believed to play a crucial role in determining the behavior of water at interfaces.39–47 Recently, important advances have been made in understanding the mechanism of water crystallization affected by confinement or electric fields.4,35,48–51 However, there is still a)Authors to whom correspondence should be addressed. Electronic ad-
dresses: [email protected]. and [email protected]. 0021-9606/2015/142(13)/134704/5/$30.00
debate about whether electric fields induce crystallization of liquid water,52,53 or electromelting of monolayer ice, and we still lack a complete understanding of the effect of electric fields and nanoconfinement on the crystallization of water. Fortunately, computer simulations can provide valuable details and information about this complex issue. In this paper, we investigate the influence of nanoscale confinement and electric fields generated by charged surfaces on the crystallization behavior of water molecules. We identified two types of 2D ordered ice: type 1 mainly induced by nanoscale confinement and type 2 mainly induced by the surface charges. The transition process between the two types of ice is controlled by the electric field generated by the charged surfaces and the diffusion coefficient of water molecules increases to the maximum just when the 2D structure of ice transits from type 1 to type 2. II. METHODOLOGY
GROMACS simulation package 4.0.754 was used to explore the phase behavior of water confined between two nanoscale surfaces (walls). The simulation framework is illustrated in Figure 1. All molecular visualizations in this paper were generated using the VMD software package.55 The model solid surfaces had the same hexagonal structure as graphite, parallel to the x- y plane. Positive and negative charges of the same magnitude 0 e ≤ q ≤ 1.0 e were assigned to atoms located diagonally in neighboring hexagons. Overall, the modeled solid surfaces were neutral. The wall–wall separation was set to 0.95 nm, which just accommodates two layers of water molecules. The space between the two walls was filled with 890 TIP4P water molecules while avoiding overlap with any existing van der Waals radii. Our simulations were performed using a time step of 2.0 fs in an NVT ensemble. Temperature (T) was maintained
142, 134704-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
134704-2
Mei et al.
J. Chem. Phys. 142, 134704 (2015)
FIG. 1. (a) Snapshot of the simulation system with bilayer ice confined between two plates with charge dipoles. (b) Geometry of the surface model. Red and blue spheres represent the atoms of walls with positive and negative charges, respectively, and the cyan spheres represent neutral atoms.
III. RESULTS AND DISCUSSION
When q increased to the range of 0.1 e ≤ q ≤ 0.3 e, the diffusion coefficient of the ice increased slowly with increasing q. Interestingly, we noted that the lateral diffusion coefficient increased sharply in the range of 0.3 e < q ≤ 0.5 e. However, it decreased dramatically when q was more than 0.5 e. As q increased further, the lateral diffusion coefficient decreased to a constant. The diffusion coefficient of q = 0.5 e was about 100 times that of q = 0 e, indicating that electromelting occurred. At q = 0.5 e, we noted that the water networks were disrupted (see the middle inset in Fig. 2), which sharply increased the mobility of water molecules. However, as q increased further, a similar but new 2D ordered hexagonal ice (see the right inset in Fig. 2), called type 2, formed. So the diffusion coefficient decreased dramatically. But at T = 300 K, the profile of the diffusion coefficient was very different from that at T = 240 K. The right top inset in Fig. 2 shows that the diffusion coefficient decreases monotonically with the charge. This is because the water molecules are almost frozen by the charge (>0.6 e) on the surfaces. Therefore, 2D ordered structures can form even at room temperature, which is consistent with our previous work.40 But without charge, ice was not found in our simulation at T = 300 K. Figures 3(a)–3(c) show the typical arrangement of water molecules confined inside the walls at T = 240 K. For the uncharged walls, the water structure was a 2D ordered hexagonal configuration resembling that of Fig. 3(a). Water molecules are connected by hydrogen bonds. There are two
To understand the effects of electric fields and nanoscale confinement on the behavior of water molecules, we calculated the lateral diffusion coefficient of the water molecules, which is an important parameter indicative of the diffusion mobility. The results are shown in Fig. 2. Our simulation results suggest that the lateral diffusion coefficient of water liquid confined inside two uncharged plates for T = 300 K was about 2.7 × 10−5 cm2 s−1 (see the inset in Fig. 2). For T = 240 K, below the bulk crystallization temperature, the diffusion coefficient of the ice was about three orders lower in magnitude than the liquid water, indicating the phase transition between liquid and solid. From the snapshot shown in the left inset in Fig. 2, 2D ordered hexagonal ice (type 1) forms, which is the typical arrangements of water molecules confined inside the walls at T = 240 K. Only a few defects were observed in our simulation results. Fig. 2 shows that the mobility of the ice remained nearly unchanged for q ≤ 0.1 e, with the lateral diffusion coefficient fluctuating slightly on the order of 10−8 cm2 s−1.
FIG. 2. In-plane diffusion coefficient of water at T = 240 K and at T = 300 K (inset) as a function of q.
using a Nose–Hoover thermostat. The carbon atoms of the walls were represented with Lennard–Jones parameters ε ss = 0.105 kcal/mol, σss = 3.343 Å. Here, we used ε ss and σss, not using ε cc and σcc, to emphasize the model system. These parameters were the same as our previous work.39,40 The simulation box size was 6.390 × 6.816 × 4.110 nm3. Periodic conditions were applied along x and y. The cutoff distance for the Lennard–Jones interactions was 15 Å. The long-range electrostatic interactions were computed by using the particle mesh Ewald method (real space cutoff, 10 Å; reciprocal space gridding, 1.2 Å, fourth-order interpolation). The surface atoms were restrained at their original positions during the simulations. To study the influence of charged surfaces on the crystallization of the confined water molecules, we simulated two layers of water molecules confined between two charged plates with d = 0.95 nm. We started the simulations from a liquid equilibrated at T = 300 K. The system was then quenched to temperature T = 240 K, below the bulk crystallization temperature (273 K) in the experiment under normal condition, and the simulation was performed until the crystallization was complete. According to the simulation studies, the TIP4P water freezes at about 214 K under normal condition,56 but it freezes at about 240 K confined inside the nanoscale space.57
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
134704-3
Mei et al.
J. Chem. Phys. 142, 134704 (2015)
FIG. 3. ((a)–(c)) Top view of typical 2D bilayer ices formed between two plates at q = 0 e, 0.5 e, 1 e. Red and white spheres represent oxygen and hydrogen atoms, and blue dotted lines represent hydrogen bonds.
types of water molecules; one type is called P-type, where the dipole of water molecule is parallel to the plates, and the other is called V-type, where the dipole of water molecule is vertical to the plates. Fig. 3(a) shows that there was a hexagonal ring consisting of six water molecules denoted by 1v, 2v, 3p, 4v, 5v, 6p in layer 1. Correspondingly, in layer 2, the hexagonal ring consisted of water molecules 1p, 2p, 3v, 4p, 5p, 6v. Interestingly, there are hydrogen bonds between P-type water molecules in the above layer (layer 1) and the V-type water molecules in the bottom layer (layer 2). As q increased, the ordered hexagonal structure was disturbed. At q = 0.5 e, the water networks were particularly disrupted (see Fig. 3(b)). However, as q increased, a similar but new 2D ordered hexagonal structure formed again. When q was sufficiently large (q > 0.7 e), half of the water molecules (P-type water molecules) bound to the charged surface mainly through attractive interactions between their O atoms and the positive charges on the surface; for the other half (V-type water molecules), any hydrogen atom that had a strong attractive interaction with a negative charge on the surface necessarily pointed towards the surface. Unlike the previous water structures, which have connecting hydrogen bonds between two water layers when q = 0 e, there were no hydrogen bonds connecting the two water layers when q = 1.0 e. The structures of the ices in layer 1 and layer 2 were the same (see Fig. 3(c)). Figure 4 shows the density profiles of the water confined inside two plates with different q. The distance between the two parallel plates was 0.95 nm and there were two layers of water molecules. Here, we plotted half of the density profile (one layer of ice) owing to symmetry. The local densities were calculated by dividing the water into many slabs in the z direction, then taking the statistical average for the density of each slab. Here, we used the coordinates of oxygen in the water molecules to calculate the density profile. From our simulation results, two sharp density peaks near the plates were observed. The location of the density peak was 0.34 nm away from the plate (shown in Fig. 4). For 0 e < q ≤ 0.3 e, the influence of the charge on the profile of the water density was negligible.
As the charge q increased further, the density peaks became flatter and the locations of the peaks moved to the charged plate. This indicated that the position of the two layers of ice changed, especially from q = 0.4 e to q = 0.5 e. At q = 0.5 e, the density profile jumped indicating that the phase transition occurred. For q > 0.5 e, the peak of water density recovered. However, the location moved even further to the charged plate. For q = 0.9 e, the location of the first ice was about 0.29 nm. Simulation results showed that the 2D ordered hexagonal structure of ice, in which water molecules are connected by strong hydrogen bonds, was controlled by the charges on the surface (see Fig. 3). Interestingly, the diffusion coefficient of water molecules increased to the maximum when the structure of ice transitioned from type 1 to type 2. The special hydrogenbonding network is responsible for the peculiar properties of water molecules confined inside the charged surfaces. In our simulation system, there were two layers (layer 1 and layer 2) of water molecules confined inside two parallel plates. We computed the average number of hydrogen bonds formed by a water molecule with its neighboring
FIG. 4. Density profile along the z axis for different q.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
134704-4
Mei et al.
water molecules in the same layer (namely, 1–1 and 2–2, respectively), as well as the average number of hydrogen bonds formed between water molecules in different layers (namely, 1–2). Here, we adopted a geometric definition of hydrogen bonds according to which water pair was hydrogen-bonded if the O–O distance was less than 3.5 Å and simultaneously the bonded O–H · · · O angle was less than 30◦.58 The average number of hydrogen bonds is shown in Fig. 5. Our simulation results showed that the average number of H-bonds between water molecules in the same layers (1–1 and 2–2) and in the different layers (1–2) remained almost constant when q < 0.3 e. However, the most noticeable change was for the number of the H-bonds between layer 1 and layer 2 (1–2), which decreased from 1.0 to 0.2 as q increased over the interval of 0.3 e – 0.7 e. The number of H-bonds in the same layers (1–1 and 2–2) changed only slightly. The total H-bonds of each water molecule decreased from 1.9 to 1.5 as q increased to 0.7 e. Note that when q was sufficiently large (q > 0.7 e), for each water molecule, there was approximately one occupied charge site to account for the electrostatic interactions from the charges on the plates. Half of the water molecules bound to the charged plates mainly through attractive Coulomb interactions between their O atoms and the positive charges on the plates, whereas the other half bound to the wall through the strong electrostatic interactions between H atoms and the negative charges on the plates (the –OH vector pointing towards the charged plates). This can also be clearly seen in Fig. 3. When q was close to 0.5 e, the average number of H-bonds between layer 1 and layer 2 decreased sharply. Correspondingly, the total average number of H-bonds formed by water with its neighboring water molecules decreased too. Meanwhile, the diffusion coefficient of the water molecules increased dramatically. The polar nature of water is one of its defining characteristics. In a water molecule, the oxygen atom and hydrogen atoms share electrons in covalent bonds, but the sharing is not equal. For the TIP4P model, q(H) is 0.52 e and q(M) is −2q(H). So, when q is close to qc = 0.5 e, the charges on the surfaces control the balance between the type 1 and type 2 ice. The probability of a water molecule binding to the charged surface is the same as the probability of the water binding to neighboring water molecules in the other layer. Consequently, the probability of finding type 1 and type 2 ice is similar. The barrier for water molecules to switch from one type to the other of ice decreases. From the above, it was concluded that the special hydrogen-bonding network balance was responsible for the abnormal increase in the diffusion coefficient of water molecules when q was close to 0.5 e. This phenomenon is attributed to the competition between two kinds of interaction; one is the interaction between water molecules, and the other is the interaction of water molecules with the charged surfaces. We know that both confinement and the special arrangement of the charge on a surface can induce the hexagonal ordered ice. The structure can be described by the angle Φ between the projection of a water dipole orientation onto the plate and the x direction. The probability distribution of Φ is shown in Fig. 6. For the neutral plates, there was a clear preference in the orientation of water dipoles with Φ = 54◦, 127◦. These peaks corresponded perfectly to the
J. Chem. Phys. 142, 134704 (2015)
FIG. 5. Average number of hydrogen bonds formed by a water molecule in different shells and between shells with respect to q.
hexagonal 2D hydrogen-bond networks. The heights of the peaks decreased as q increased. At q = 0.5 e, the heights of the peaks were at their lowest and the 2D ordered hexagonal structure was disturbed. The average number of hydrogen bonds decreased to the minimum, which increased the diffusion coefficient. However, the profile of the probability distribution of Φ did not become flat, indicating that the structure of water molecules was still partly ordered. As q > 0.5 e, the peaks increased in amplitude with q. The new 2D ordered structure (type 2) was formed and is consistent with the conclusions from analyzing the simulation results shown in Figs. 3 and 5. For our system, the electric field generated by the charged surfaces is very different from the uniform electric. The magnitude and direction of the electric field vary with the position due to the hexagonal pattern of the charged walls. It is not easy to rearrange the structure of water molecules confined inside the nanoscale space only by introducing the external uniform electronic field due to strong hydrogen bonds between water molecules. It is known that the asymmetry of the water molecule leads to a dipole moment, which is about 2.18 D for TIP4P model.59 The interaction energy between
FIG. 6. Probability distribution of the angle Φ, where Φ is the angle formed between the x–y plane projection of a water molecule dipole orientation and the x direction.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
134704-5
Mei et al.
two water molecules is about 6.24 kcal/mol. From Fig. 5, we can see that one hydrogen bond will be broken when the 2D ordered ice transforms from type 1 to type 2. So we can calculate the magnitude of the critical electric field in which water structure is just rearranged (polarized) according to equation E = U/2p, where E is electric field, U is the interaction energy between two water molecules, and p is the dipole moment of a water molecule. So, in this case of TIP4P water model, the critical electric field is about 2.94 × 109 V/m. For our simulation system, the critical charge value is about 0.5 e. And the distance between water molecules and the wall ´ Then, the magnitude of the electronic field is about 2.9 Å. in the position of the water molecule above the charge is about 2.14 × 109 V/m. The estimated value is consistent with the result in Ref. 35. It is still difficult to produce so huge electronic field strength (109 V/m) in experiment. However, we can introduce charges to control the dipoles’ directions of the water molecules. The corresponding realistic system could be Talc, InSb(110) surface, hydroxylated alumina, and hydroxylated silica surfaces. Beside the corresponding realistic system, there are many experimental methods of modifying the walls with charged-group in the appropriate positions, such as intentional doping the functional groups of –COO− and –NH3+. In fact, the similar mechanism can be found in the antifreeze proteins.
IV. CONCLUSIONS
In this work, we systematically investigated the effects of charged surfaces on the phase behavior of two layers of water molecules confined inside a nanoscale space. We found that charges on surfaces greatly influence the configuration and dynamics behavior of ice layers, specifically the structure and the stability of the two-dimensional hexagonal hydrogen-bond network. The in-plane diffusion coefficient of water molecules increased dramatically when q was close to qc = 0.5 e due to the transformation of the hydrogen bond network. The mechanism behind this unexpected phenomenon rests in the competition between two kinds of interaction; one is the interaction between water molecules, and the other is the interaction of water molecules with the charged surfaces.
ACKNOWLEDGMENTS
We would like to thank Professor Haiping Fang for helpful discussions. This work is partially supported by the National Natural Science Foundation of China (Grant Nos. 11405146, 11204341, and 61274099), and Zhejiang Provincial Science and Technology Key Innovation Team No. 2011R50012 and Zhejiang Provincial Key Laboratory No. 2013E10022, the China Scholarship Council (No. 201308330031), the Knowledge Innovation Program of SINAP, the Knowledge Innovation Program of the Chinese Academy of Sciences, Shanghai Supercomputer Center of China, Deepcomp 7000 and ScGrid
J. Chem. Phys. 142, 134704 (2015)
of Supercomputing Center, and Computer Network Information Center of Chinese Academy of Sciences. 1J.
P. Doye et al., Phys. Chem. Chem. Phys. 9, 2197 (2007). Lv et al., ACS Nano 8, 3152 (2014). 3M. J. Kuiper, P. L. Davies, and V. K. Walker, Biophys. J. 81, 3560 (2001). 4G. Sun et al., Soft Matter 9, 11374 (2013). 5H. Kiani and D.-W. Sun, Trends Food Sci. Technol. 22, 407 (2011). 6K. Meister et al., Proc. Natl. Acad. Sci. U. S. A. 110, 1617 (2013). 7N. Kastelowitz, J. C. Johnston, and V. Molinero, J. Chem. Phys. 132, 124511 (2010). 8J. Bai and X. C. Zeng, Proc. Natl. Acad. Sci. U. S. A. 109, 21240 (2012). 9M. Garcia-Ratés et al., J. Phys. Chem. B 115, 5980 (2011). 10R. Srivastava et al., J. Phys. Chem. C 115, 12448 (2011). 11S. Han et al., Nature Phys. 6, 685 (2010). 12N. Giovambattista, P. J. Rossky, and P. G. Debenedetti, Phys. Rev. Lett. 102, 50603 (2009). 13D. Takaiwa et al., Proc. Natl. Acad. Sci. U. S. A. 105, 39 (2008). 14J. M. Zanotti, M. C. Bellissent-Funel, and A. I. Kolesnikov, Eur. Phys. J.: Spec. Top. 141, 227 (2007). 15K. Koga and H. Tanaka, J. Chem. Phys. 122, 104711 (2005). 16K. Koga, H. Tanaka, and X. Zeng, Nature 408, 564 (2000). 17E. G. Solveyra et al., J. Phys. Chem. B 115, 14196 (2011). 18B. K. Agrawal et al., Phys. Rev. B 75, 195420 (2007). 19P. Levitz et al., Soft Matter 9, 8654 (2013). 20Y. Long et al., Phys. Chem. Chem. Phys. 13, 17163 (2011). 21F. G. Alabarse et al., Phys. Rev. Lett. 109, 035701 (2012). 22M. Odelius, M. Bernasconi, and M. Parrinello, Phys. Rev. Lett. 78, 2855 (1997). 23P. B. Miranda et al., Phys. Rev. Lett. 81, 5876 (1998). 24R. J. Mashl et al., Nano Lett. 3, 589 (2003). 25L. B. Krott and M. C. Barbosa, J. Chem. Phys. 138, 084505 (2013). 26W. Zhao et al., Chem. Sci. 5, 1757 (2014). 27K. Koga et al., Nature 412, 802 (2001). 28K. Koga et al., J. Chem. Phys. 113, 5037 (2000). 29D. Takaiwa, K. Koga, and H. Tanaka, Mol. Simul. 33, 127 (2007). 30J. Shiomi, T. Kimura, and S. Maruyama, J. Phys. Chem. C 111, 12188 (2007). 31T. Kurita, S. Okada, and A. Oshiyama, Phys. Rev. B 75, 205424 (2007). 32H.-J. Lu et al., Chin. Phys. Lett. 25, 1145 (2008). 33C. T. Kresge et al., Nature 359, 710 (1992). 34M. Erko et al., Phys. Rev. B 84, 104205 (2011). 35H. Qiu and W. L. Guo, Phys. Rev. Lett. 110, 195701 (2013). 36R. Zangi and A. E. Mark, Phys. Rev. Lett. 91, 025502 (2003). 37R. Zangi and A. E. Mark, J. Chem. Phys. 119, 1694 (2003). 38J. Slovak et al., Phys. Rev. E 60, 5833 (1999). 39X. Zhou et al., J. Chem. Phys. 138, 204710 (2013). 40C. Wang et al., Phys. Rev. Lett. 103, 137801 (2009). 41N. Giovambattista, P. G. Debenedetti, and P. J. Rossky, J. Phys. Chem. B 111, 9581 (2007). 42S. R.-V. Castrillón et al., J. Phys. Chem. B 113, 1438 (2009). 43K. Bagchi and S. Roy, J. Phys. Chem. B 118, 3805 (2014). 44C. Wang, Y. Yang, and H. Fang, Sci. China: Phys., Mech. Astron. 57, 802 (2014). 45C. L. Wang, J. Y. Li, and H. P. Fang, Rendi. Lincei 22, 1 (2011). 46C. L. Wang et al., Sci. Rep. 2, 358 (2012). 47C. L. Wang et al., J. Phys. Chem. C 115, 3018 (2011). 48Z. Fu et al., J. Chem. Phys. 134, 154507 (2011). 49P. Kumar and S. Han, J. Chem. Phys. 137, 114510 (2012). 50Z. Qian and G. Wei, J. Phys. Chem. A 118, 8922 (2014). 51X. Zhu, Q. Yuan, and Y.-P. Zhao, Nanoscale 6, 5432 (2014). 52I. M. Svishchev and P. G. Kusalik, J. Am. Chem. Soc. 118, 649 (1996). 53M. Gavish et al., Science 256, 815 (1992). 54B. Hess et al., J. Chem. Theory Comput. 4, 435 (2008). 55W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graphics 14, 33 (1996). 56M. J. Vlot, J. Huinink, and J. P. van der Eerden, J. Chem. Phys. 110, 55 (1999). 57K. Koga, X. C. Zeng, and H. Tanaka, Phys. Rev. Lett. 79, 5262 (1997). 58H. F. Xu and B. J. Berne, J. Phys. Chem. B 105, 11929 (2001). 59W. L. Jorgensen et al., J. Chem. Phys. 79, 926 (1983). 2J.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.
A transition between bistable ice when coupling electric field and nanoconfinement Feng Mei,1,2 Xiaoyan Zhou,1,3 Jianlong Kou,1 Fengmin Wu,1,3 Chunlei Wang,4,a) and Hangjun Lu1,a) 1
Department of Physics, Zhejiang Normal University, Jinhua 321004, China Department of Physics, Zhejiang University, Hangzhou 310027, China 3 Department of Physics and Institute of Theoretical Physics, Shanxi University, Taiyuan 030006, China 4 Division of Interfacial Water and Key Laboratory of Interfacial Physics and Technology, Shanghai Institute of Applied Physics, Chinese Academy of Sciences, P.O. Box 800-204, Shanghai 201800, China 2
(Received 16 January 2015; accepted 17 March 2015; published online 2 April 2015) The effects of an electric field on the phase behavior of water confined inside a nanoscale space were studied using molecular dynamics simulations. It was found that the diffusion coefficient of water reaches its maximum when value of the surfaces’ charge is at the threshold, qc = 0.5e. This unexpected phenomenon was attributed to the intermediate state between two stable ice states induced by nanoconfinement and the electric field generated by charged surfaces, respectively. Our finding is helpful to understand electromelting and electrofreezing of water under nanoconfinement with the electric field. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4916521]
I. INTRODUCTION
Understanding, controlling, and predicting the crystallization process of nano-confined liquids are fundamental to physics, chemistry, biology, and engineering applications.1–5 Controlling crystallization in cells can determine the ability of many organisms such as polar fish, insects, bacteria, and plants, to survive and thrive in subfreezing habitats.6 In the food, pharmaceutical and chemical industries, controlling crystallization in confined spaces is critical for both functional performance and regulatory compliance.4 However, the crystallization behavior of water is a complex phase transition that presents a challenge to statistical mechanics. Confinement of water inside nanoscale geometries or in the vicinity of surfaces can lead to interesting properties and increase the richness of its phase behavior, which differs considerably from bulk systems.7–26 A number of additional ice phases have been observed by simulating water confined within carbon nanotubes.7,13,16,18,27–32 Experimental results showed that decreasing the pore diameter depressed the freezing/melting point according to a modified Gibbs–Thomson equation.33,34 Computational and experimental studies have demonstrated the spontaneous freezing of monolayer, bilayer, and three-layer ice confined between parallel plates separated by nanometer gaps.7,26,35–38 Interestingly, there is still a debate about whether the nanoscale confinement can favor or suppress crystallization. Another factor that may affect water crystallization is surface polarity, which is believed to play a crucial role in determining the behavior of water at interfaces.39–47 Recently, important advances have been made in understanding the mechanism of water crystallization affected by confinement or electric fields.4,35,48–51 However, there is still a)Authors to whom correspondence should be addressed. Electronic ad-
dresses: [email protected]. and [email protected]. 0021-9606/2015/142(13)/134704/5/$30.00
debate about whether electric fields induce crystallization of liquid water,52,53 or electromelting of monolayer ice, and we still lack a complete understanding of the effect of electric fields and nanoconfinement on the crystallization of water. Fortunately, computer simulations can provide valuable details and information about this complex issue. In this paper, we investigate the influence of nanoscale confinement and electric fields generated by charged surfaces on the crystallization behavior of water molecules. We identified two types of 2D ordered ice: type 1 mainly induced by nanoscale confinement and type 2 mainly induced by the surface charges. The transition process between the two types of ice is controlled by the electric field generated by the charged surfaces and the diffusion coefficient of water molecules increases to the maximum just when the 2D structure of ice transits from type 1 to type 2. II. METHODOLOGY
GROMACS simulation package 4.0.754 was used to explore the phase behavior of water confined between two nanoscale surfaces (walls). The simulation framework is illustrated in Figure 1. All molecular visualizations in this paper were generated using the VMD software package.55 The model solid surfaces had the same hexagonal structure as graphite, parallel to the x- y plane. Positive and negative charges of the same magnitude 0 e ≤ q ≤ 1.0 e were assigned to atoms located diagonally in neighboring hexagons. Overall, the modeled solid surfaces were neutral. The wall–wall separation was set to 0.95 nm, which just accommodates two layers of water molecules. The space between the two walls was filled with 890 TIP4P water molecules while avoiding overlap with any existing van der Waals radii. Our simulations were performed using a time step of 2.0 fs in an NVT ensemble. Temperature (T) was maintained
142, 134704-1
© 2015 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
134704-2
Mei et al.
J. Chem. Phys. 142, 134704 (2015)
FIG. 1. (a) Snapshot of the simulation system with bilayer ice confined between two plates with charge dipoles. (b) Geometry of the surface model. Red and blue spheres represent the atoms of walls with positive and negative charges, respectively, and the cyan spheres represent neutral atoms.
III. RESULTS AND DISCUSSION
When q increased to the range of 0.1 e ≤ q ≤ 0.3 e, the diffusion coefficient of the ice increased slowly with increasing q. Interestingly, we noted that the lateral diffusion coefficient increased sharply in the range of 0.3 e < q ≤ 0.5 e. However, it decreased dramatically when q was more than 0.5 e. As q increased further, the lateral diffusion coefficient decreased to a constant. The diffusion coefficient of q = 0.5 e was about 100 times that of q = 0 e, indicating that electromelting occurred. At q = 0.5 e, we noted that the water networks were disrupted (see the middle inset in Fig. 2), which sharply increased the mobility of water molecules. However, as q increased further, a similar but new 2D ordered hexagonal ice (see the right inset in Fig. 2), called type 2, formed. So the diffusion coefficient decreased dramatically. But at T = 300 K, the profile of the diffusion coefficient was very different from that at T = 240 K. The right top inset in Fig. 2 shows that the diffusion coefficient decreases monotonically with the charge. This is because the water molecules are almost frozen by the charge (>0.6 e) on the surfaces. Therefore, 2D ordered structures can form even at room temperature, which is consistent with our previous work.40 But without charge, ice was not found in our simulation at T = 300 K. Figures 3(a)–3(c) show the typical arrangement of water molecules confined inside the walls at T = 240 K. For the uncharged walls, the water structure was a 2D ordered hexagonal configuration resembling that of Fig. 3(a). Water molecules are connected by hydrogen bonds. There are two
To understand the effects of electric fields and nanoscale confinement on the behavior of water molecules, we calculated the lateral diffusion coefficient of the water molecules, which is an important parameter indicative of the diffusion mobility. The results are shown in Fig. 2. Our simulation results suggest that the lateral diffusion coefficient of water liquid confined inside two uncharged plates for T = 300 K was about 2.7 × 10−5 cm2 s−1 (see the inset in Fig. 2). For T = 240 K, below the bulk crystallization temperature, the diffusion coefficient of the ice was about three orders lower in magnitude than the liquid water, indicating the phase transition between liquid and solid. From the snapshot shown in the left inset in Fig. 2, 2D ordered hexagonal ice (type 1) forms, which is the typical arrangements of water molecules confined inside the walls at T = 240 K. Only a few defects were observed in our simulation results. Fig. 2 shows that the mobility of the ice remained nearly unchanged for q ≤ 0.1 e, with the lateral diffusion coefficient fluctuating slightly on the order of 10−8 cm2 s−1.
FIG. 2. In-plane diffusion coefficient of water at T = 240 K and at T = 300 K (inset) as a function of q.
using a Nose–Hoover thermostat. The carbon atoms of the walls were represented with Lennard–Jones parameters ε ss = 0.105 kcal/mol, σss = 3.343 Å. Here, we used ε ss and σss, not using ε cc and σcc, to emphasize the model system. These parameters were the same as our previous work.39,40 The simulation box size was 6.390 × 6.816 × 4.110 nm3. Periodic conditions were applied along x and y. The cutoff distance for the Lennard–Jones interactions was 15 Å. The long-range electrostatic interactions were computed by using the particle mesh Ewald method (real space cutoff, 10 Å; reciprocal space gridding, 1.2 Å, fourth-order interpolation). The surface atoms were restrained at their original positions during the simulations. To study the influence of charged surfaces on the crystallization of the confined water molecules, we simulated two layers of water molecules confined between two charged plates with d = 0.95 nm. We started the simulations from a liquid equilibrated at T = 300 K. The system was then quenched to temperature T = 240 K, below the bulk crystallization temperature (273 K) in the experiment under normal condition, and the simulation was performed until the crystallization was complete. According to the simulation studies, the TIP4P water freezes at about 214 K under normal condition,56 but it freezes at about 240 K confined inside the nanoscale space.57
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
134704-3
Mei et al.
J. Chem. Phys. 142, 134704 (2015)
FIG. 3. ((a)–(c)) Top view of typical 2D bilayer ices formed between two plates at q = 0 e, 0.5 e, 1 e. Red and white spheres represent oxygen and hydrogen atoms, and blue dotted lines represent hydrogen bonds.
types of water molecules; one type is called P-type, where the dipole of water molecule is parallel to the plates, and the other is called V-type, where the dipole of water molecule is vertical to the plates. Fig. 3(a) shows that there was a hexagonal ring consisting of six water molecules denoted by 1v, 2v, 3p, 4v, 5v, 6p in layer 1. Correspondingly, in layer 2, the hexagonal ring consisted of water molecules 1p, 2p, 3v, 4p, 5p, 6v. Interestingly, there are hydrogen bonds between P-type water molecules in the above layer (layer 1) and the V-type water molecules in the bottom layer (layer 2). As q increased, the ordered hexagonal structure was disturbed. At q = 0.5 e, the water networks were particularly disrupted (see Fig. 3(b)). However, as q increased, a similar but new 2D ordered hexagonal structure formed again. When q was sufficiently large (q > 0.7 e), half of the water molecules (P-type water molecules) bound to the charged surface mainly through attractive interactions between their O atoms and the positive charges on the surface; for the other half (V-type water molecules), any hydrogen atom that had a strong attractive interaction with a negative charge on the surface necessarily pointed towards the surface. Unlike the previous water structures, which have connecting hydrogen bonds between two water layers when q = 0 e, there were no hydrogen bonds connecting the two water layers when q = 1.0 e. The structures of the ices in layer 1 and layer 2 were the same (see Fig. 3(c)). Figure 4 shows the density profiles of the water confined inside two plates with different q. The distance between the two parallel plates was 0.95 nm and there were two layers of water molecules. Here, we plotted half of the density profile (one layer of ice) owing to symmetry. The local densities were calculated by dividing the water into many slabs in the z direction, then taking the statistical average for the density of each slab. Here, we used the coordinates of oxygen in the water molecules to calculate the density profile. From our simulation results, two sharp density peaks near the plates were observed. The location of the density peak was 0.34 nm away from the plate (shown in Fig. 4). For 0 e < q ≤ 0.3 e, the influence of the charge on the profile of the water density was negligible.
As the charge q increased further, the density peaks became flatter and the locations of the peaks moved to the charged plate. This indicated that the position of the two layers of ice changed, especially from q = 0.4 e to q = 0.5 e. At q = 0.5 e, the density profile jumped indicating that the phase transition occurred. For q > 0.5 e, the peak of water density recovered. However, the location moved even further to the charged plate. For q = 0.9 e, the location of the first ice was about 0.29 nm. Simulation results showed that the 2D ordered hexagonal structure of ice, in which water molecules are connected by strong hydrogen bonds, was controlled by the charges on the surface (see Fig. 3). Interestingly, the diffusion coefficient of water molecules increased to the maximum when the structure of ice transitioned from type 1 to type 2. The special hydrogenbonding network is responsible for the peculiar properties of water molecules confined inside the charged surfaces. In our simulation system, there were two layers (layer 1 and layer 2) of water molecules confined inside two parallel plates. We computed the average number of hydrogen bonds formed by a water molecule with its neighboring
FIG. 4. Density profile along the z axis for different q.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
134704-4
Mei et al.
water molecules in the same layer (namely, 1–1 and 2–2, respectively), as well as the average number of hydrogen bonds formed between water molecules in different layers (namely, 1–2). Here, we adopted a geometric definition of hydrogen bonds according to which water pair was hydrogen-bonded if the O–O distance was less than 3.5 Å and simultaneously the bonded O–H · · · O angle was less than 30◦.58 The average number of hydrogen bonds is shown in Fig. 5. Our simulation results showed that the average number of H-bonds between water molecules in the same layers (1–1 and 2–2) and in the different layers (1–2) remained almost constant when q < 0.3 e. However, the most noticeable change was for the number of the H-bonds between layer 1 and layer 2 (1–2), which decreased from 1.0 to 0.2 as q increased over the interval of 0.3 e – 0.7 e. The number of H-bonds in the same layers (1–1 and 2–2) changed only slightly. The total H-bonds of each water molecule decreased from 1.9 to 1.5 as q increased to 0.7 e. Note that when q was sufficiently large (q > 0.7 e), for each water molecule, there was approximately one occupied charge site to account for the electrostatic interactions from the charges on the plates. Half of the water molecules bound to the charged plates mainly through attractive Coulomb interactions between their O atoms and the positive charges on the plates, whereas the other half bound to the wall through the strong electrostatic interactions between H atoms and the negative charges on the plates (the –OH vector pointing towards the charged plates). This can also be clearly seen in Fig. 3. When q was close to 0.5 e, the average number of H-bonds between layer 1 and layer 2 decreased sharply. Correspondingly, the total average number of H-bonds formed by water with its neighboring water molecules decreased too. Meanwhile, the diffusion coefficient of the water molecules increased dramatically. The polar nature of water is one of its defining characteristics. In a water molecule, the oxygen atom and hydrogen atoms share electrons in covalent bonds, but the sharing is not equal. For the TIP4P model, q(H) is 0.52 e and q(M) is −2q(H). So, when q is close to qc = 0.5 e, the charges on the surfaces control the balance between the type 1 and type 2 ice. The probability of a water molecule binding to the charged surface is the same as the probability of the water binding to neighboring water molecules in the other layer. Consequently, the probability of finding type 1 and type 2 ice is similar. The barrier for water molecules to switch from one type to the other of ice decreases. From the above, it was concluded that the special hydrogen-bonding network balance was responsible for the abnormal increase in the diffusion coefficient of water molecules when q was close to 0.5 e. This phenomenon is attributed to the competition between two kinds of interaction; one is the interaction between water molecules, and the other is the interaction of water molecules with the charged surfaces. We know that both confinement and the special arrangement of the charge on a surface can induce the hexagonal ordered ice. The structure can be described by the angle Φ between the projection of a water dipole orientation onto the plate and the x direction. The probability distribution of Φ is shown in Fig. 6. For the neutral plates, there was a clear preference in the orientation of water dipoles with Φ = 54◦, 127◦. These peaks corresponded perfectly to the
J. Chem. Phys. 142, 134704 (2015)
FIG. 5. Average number of hydrogen bonds formed by a water molecule in different shells and between shells with respect to q.
hexagonal 2D hydrogen-bond networks. The heights of the peaks decreased as q increased. At q = 0.5 e, the heights of the peaks were at their lowest and the 2D ordered hexagonal structure was disturbed. The average number of hydrogen bonds decreased to the minimum, which increased the diffusion coefficient. However, the profile of the probability distribution of Φ did not become flat, indicating that the structure of water molecules was still partly ordered. As q > 0.5 e, the peaks increased in amplitude with q. The new 2D ordered structure (type 2) was formed and is consistent with the conclusions from analyzing the simulation results shown in Figs. 3 and 5. For our system, the electric field generated by the charged surfaces is very different from the uniform electric. The magnitude and direction of the electric field vary with the position due to the hexagonal pattern of the charged walls. It is not easy to rearrange the structure of water molecules confined inside the nanoscale space only by introducing the external uniform electronic field due to strong hydrogen bonds between water molecules. It is known that the asymmetry of the water molecule leads to a dipole moment, which is about 2.18 D for TIP4P model.59 The interaction energy between
FIG. 6. Probability distribution of the angle Φ, where Φ is the angle formed between the x–y plane projection of a water molecule dipole orientation and the x direction.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
134704-5
Mei et al.
two water molecules is about 6.24 kcal/mol. From Fig. 5, we can see that one hydrogen bond will be broken when the 2D ordered ice transforms from type 1 to type 2. So we can calculate the magnitude of the critical electric field in which water structure is just rearranged (polarized) according to equation E = U/2p, where E is electric field, U is the interaction energy between two water molecules, and p is the dipole moment of a water molecule. So, in this case of TIP4P water model, the critical electric field is about 2.94 × 109 V/m. For our simulation system, the critical charge value is about 0.5 e. And the distance between water molecules and the wall ´ Then, the magnitude of the electronic field is about 2.9 Å. in the position of the water molecule above the charge is about 2.14 × 109 V/m. The estimated value is consistent with the result in Ref. 35. It is still difficult to produce so huge electronic field strength (109 V/m) in experiment. However, we can introduce charges to control the dipoles’ directions of the water molecules. The corresponding realistic system could be Talc, InSb(110) surface, hydroxylated alumina, and hydroxylated silica surfaces. Beside the corresponding realistic system, there are many experimental methods of modifying the walls with charged-group in the appropriate positions, such as intentional doping the functional groups of –COO− and –NH3+. In fact, the similar mechanism can be found in the antifreeze proteins.
IV. CONCLUSIONS
In this work, we systematically investigated the effects of charged surfaces on the phase behavior of two layers of water molecules confined inside a nanoscale space. We found that charges on surfaces greatly influence the configuration and dynamics behavior of ice layers, specifically the structure and the stability of the two-dimensional hexagonal hydrogen-bond network. The in-plane diffusion coefficient of water molecules increased dramatically when q was close to qc = 0.5 e due to the transformation of the hydrogen bond network. The mechanism behind this unexpected phenomenon rests in the competition between two kinds of interaction; one is the interaction between water molecules, and the other is the interaction of water molecules with the charged surfaces.
ACKNOWLEDGMENTS
We would like to thank Professor Haiping Fang for helpful discussions. This work is partially supported by the National Natural Science Foundation of China (Grant Nos. 11405146, 11204341, and 61274099), and Zhejiang Provincial Science and Technology Key Innovation Team No. 2011R50012 and Zhejiang Provincial Key Laboratory No. 2013E10022, the China Scholarship Council (No. 201308330031), the Knowledge Innovation Program of SINAP, the Knowledge Innovation Program of the Chinese Academy of Sciences, Shanghai Supercomputer Center of China, Deepcomp 7000 and ScGrid
J. Chem. Phys. 142, 134704 (2015)
of Supercomputing Center, and Computer Network Information Center of Chinese Academy of Sciences. 1J.
P. Doye et al., Phys. Chem. Chem. Phys. 9, 2197 (2007). Lv et al., ACS Nano 8, 3152 (2014). 3M. J. Kuiper, P. L. Davies, and V. K. Walker, Biophys. J. 81, 3560 (2001). 4G. Sun et al., Soft Matter 9, 11374 (2013). 5H. Kiani and D.-W. Sun, Trends Food Sci. Technol. 22, 407 (2011). 6K. Meister et al., Proc. Natl. Acad. Sci. U. S. A. 110, 1617 (2013). 7N. Kastelowitz, J. C. Johnston, and V. Molinero, J. Chem. Phys. 132, 124511 (2010). 8J. Bai and X. C. Zeng, Proc. Natl. Acad. Sci. U. S. A. 109, 21240 (2012). 9M. Garcia-Ratés et al., J. Phys. Chem. B 115, 5980 (2011). 10R. Srivastava et al., J. Phys. Chem. C 115, 12448 (2011). 11S. Han et al., Nature Phys. 6, 685 (2010). 12N. Giovambattista, P. J. Rossky, and P. G. Debenedetti, Phys. Rev. Lett. 102, 50603 (2009). 13D. Takaiwa et al., Proc. Natl. Acad. Sci. U. S. A. 105, 39 (2008). 14J. M. Zanotti, M. C. Bellissent-Funel, and A. I. Kolesnikov, Eur. Phys. J.: Spec. Top. 141, 227 (2007). 15K. Koga and H. Tanaka, J. Chem. Phys. 122, 104711 (2005). 16K. Koga, H. Tanaka, and X. Zeng, Nature 408, 564 (2000). 17E. G. Solveyra et al., J. Phys. Chem. B 115, 14196 (2011). 18B. K. Agrawal et al., Phys. Rev. B 75, 195420 (2007). 19P. Levitz et al., Soft Matter 9, 8654 (2013). 20Y. Long et al., Phys. Chem. Chem. Phys. 13, 17163 (2011). 21F. G. Alabarse et al., Phys. Rev. Lett. 109, 035701 (2012). 22M. Odelius, M. Bernasconi, and M. Parrinello, Phys. Rev. Lett. 78, 2855 (1997). 23P. B. Miranda et al., Phys. Rev. Lett. 81, 5876 (1998). 24R. J. Mashl et al., Nano Lett. 3, 589 (2003). 25L. B. Krott and M. C. Barbosa, J. Chem. Phys. 138, 084505 (2013). 26W. Zhao et al., Chem. Sci. 5, 1757 (2014). 27K. Koga et al., Nature 412, 802 (2001). 28K. Koga et al., J. Chem. Phys. 113, 5037 (2000). 29D. Takaiwa, K. Koga, and H. Tanaka, Mol. Simul. 33, 127 (2007). 30J. Shiomi, T. Kimura, and S. Maruyama, J. Phys. Chem. C 111, 12188 (2007). 31T. Kurita, S. Okada, and A. Oshiyama, Phys. Rev. B 75, 205424 (2007). 32H.-J. Lu et al., Chin. Phys. Lett. 25, 1145 (2008). 33C. T. Kresge et al., Nature 359, 710 (1992). 34M. Erko et al., Phys. Rev. B 84, 104205 (2011). 35H. Qiu and W. L. Guo, Phys. Rev. Lett. 110, 195701 (2013). 36R. Zangi and A. E. Mark, Phys. Rev. Lett. 91, 025502 (2003). 37R. Zangi and A. E. Mark, J. Chem. Phys. 119, 1694 (2003). 38J. Slovak et al., Phys. Rev. E 60, 5833 (1999). 39X. Zhou et al., J. Chem. Phys. 138, 204710 (2013). 40C. Wang et al., Phys. Rev. Lett. 103, 137801 (2009). 41N. Giovambattista, P. G. Debenedetti, and P. J. Rossky, J. Phys. Chem. B 111, 9581 (2007). 42S. R.-V. Castrillón et al., J. Phys. Chem. B 113, 1438 (2009). 43K. Bagchi and S. Roy, J. Phys. Chem. B 118, 3805 (2014). 44C. Wang, Y. Yang, and H. Fang, Sci. China: Phys., Mech. Astron. 57, 802 (2014). 45C. L. Wang, J. Y. Li, and H. P. Fang, Rendi. Lincei 22, 1 (2011). 46C. L. Wang et al., Sci. Rep. 2, 358 (2012). 47C. L. Wang et al., J. Phys. Chem. C 115, 3018 (2011). 48Z. Fu et al., J. Chem. Phys. 134, 154507 (2011). 49P. Kumar and S. Han, J. Chem. Phys. 137, 114510 (2012). 50Z. Qian and G. Wei, J. Phys. Chem. A 118, 8922 (2014). 51X. Zhu, Q. Yuan, and Y.-P. Zhao, Nanoscale 6, 5432 (2014). 52I. M. Svishchev and P. G. Kusalik, J. Am. Chem. Soc. 118, 649 (1996). 53M. Gavish et al., Science 256, 815 (1992). 54B. Hess et al., J. Chem. Theory Comput. 4, 435 (2008). 55W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graphics 14, 33 (1996). 56M. J. Vlot, J. Huinink, and J. P. van der Eerden, J. Chem. Phys. 110, 55 (1999). 57K. Koga, X. C. Zeng, and H. Tanaka, Phys. Rev. Lett. 79, 5262 (1997). 58H. F. Xu and B. J. Berne, J. Phys. Chem. B 105, 11929 (2001). 59W. L. Jorgensen et al., J. Chem. Phys. 79, 926 (1983). 2J.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 202.177.173.189 On: Wed, 08 Apr 2015 10:14:42
Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.
