The Ice−Vapor Interface and the Melting Point of Ice Ih for the

ARTICLE pubs.acs.org/JPCA

The Ice
Vapor Interface and the Melting Point of Ice Ih for the Polarizable POL3 Water Model Eva Muchov
a,† Ivan Gladich,‡ Sylvain Picaud,§ Paul N. M. Hoang,§ and Martina Roeselov
a*,‡ †

Faculty of Science, Charles University Prague, Albertov 6, Prague, 128 43, Czech Republic Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, Flemingovo n
am. 2, Prague, 166 10, Czech Republic § Institut UTINAM -UMR CNRS 6213, Universit
e de Franche-Comt
e, 25030 Besanc- on Cedex, France ‡

ABSTRACT: We use molecular dynamics simulations to determine the melting point of ice Ih for the polarizable POL3 water force field (Dang, L. X. J. Chem. Phys. 1992, 97, 2659). Simulations are performed on a slab of ice Ih with two free surfaces at several different temperatures. The analysis of the time evolution of the total energy in the course of the simulations at the set of temperatures yields the melting point of the POL3 model to be Tm = 180 ( 10 K. Moreover, the results of the simulations show that the degree of hydrogen-bond disorder occurring in the bulk of POL3 ice is larger (at the corresponding degree of undercooling) than in ice modeled by nonpolarizable water models. These results demonstrate that the POL3 water force field is rather a poor model for studying ice and ice
liquid or ice
vapor interfaces. While a number of polarizable water models have been developed over the past years, little is known about their performance in simulations of supercooled water and ice. This study thus highlights the need for testing of the existing polarizable water models over a broad range of temperatures, pressures, and phases, and developing a new polarizable water force field, reliable over larger areas of the phase diagram.

1. INTRODUCTION Over the last four decades since the first Monte Carlo and molecular dynamics studies of liquid water,1,2 computer simulations based on empirical models (force fields) have proven very useful in understanding the properties of water in molecular detail. The vast majority of computational studies have focused on liquid water under ambient conditions, exploring the properties of both water as a pure liquid and water as a solvent. However, in many cases of interest, for example, in geosciences, atmospheric sciences, as well as in engineering applications, water occurs far from room temperature and pressure. There has thus been a growing interest in employing computer simulations to gain molecular level information on water over a wide range of temperatures, pressures, and phases. In particular, considerable attention has recently been given to water in its solid form as ice or, more generally, to understanding the behavior of water at low temperatures—supercooled water and the process of ice nucleation and freezing, amorphous solid water, transitions between the numerous solid phases, etc. A particularly active area of research motivated largely by the needs of the atmospheric science community is currently focused on the properties of the ice
liquid and ice
vapor interfaces, both pure and when impurities and adsorbates are present. Indeed, the upper troposphere is characterized by low temperatures ranging from about 188 to 233 K, allowing the formation of cirrus clouds covering a large portion (up to 25%) of the Earth’s surface.3 Below 233 K, these cirrus clouds are mainly made of ice r 2011 American Chemical Society

particles,3,4 whereas, above this temperature, supercooled water droplets have also been observed.5 The interaction between trace gases and these ice particles or supercooled droplets is strongly suspected to promote heterogeneous chemistry, and therefore, they can offer a favorable environment for fast physical and chemical transformations that can modify both trace gas and aerosol compositions in the troposphere.6 In addition, the presence of ionic impurities or organic pollutants may substantially alter the properties of ice particles and supercooled droplets in cirrus clouds;7,8 however, the corresponding processes are presently not well established.7
9 It has been shown in the past decade that numerical simulations are a powerful tool for characterizing the adsorption of small organic compounds at the surface of pure ice and theoretical studies based either on the molecular dynamics technics10
13 or on the Monte Carlo methods14
19 have shown very good agreement with experimental data recorded under similar conditions. However, as far as we know, no similar studies have been performed on ice doped with ionic impurities such as Naþ, Cl
, Br
, or NO3
, to name just a few of the atmospherically relevant ones. Nevertheless, given the complexities of the naturally occurring substrates, we Special Issue: Victoria Buch Memorial Received: October 30, 2010 Revised: March 18, 2011 Published: March 31, 2011 5973

dx.doi.org/10.1021/jp110391q | J. Phys. Chem. A 2011, 115, 5973–5982

The Journal of Physical Chemistry A

ARTICLE

Table 1. Melting Point Tm (K) of Selected Water Models

a

SPC

SPC/E

TIP3P

TIP4P

TIP4P-Ew

TIP4P/2005

190.5a

215b

146a

230b

243.5b

250b

TIP4P/Ice

TIP5P

NE6

TIP4P-FQc

270b

272.5d

289d

303e

TTM3-Fc

POL3c

250f

180g

36 b

Tm obtained from free energy calculations. Tm obtained as an average of results from free energy calculations, calculations based on coexistence of the solid
liquid interface, and simulations of ice with a free surface.38 c Polarizable force field. d Tm obtained as an average of results from free energy calculations and calculations based on coexistence of the solid
liquid interface.34,35 e Reference 64. f Tm obtained from classical simulation.66 g This work. Tm obtained from simulations of ice with a free surface.

cannot hope to obtain at least a qualitatively correct understanding of important heterogeneous atmospheric chemistry and (micro)physics without taking into account ionic constituents on ice surfaces. From the point of view of classical molecular simulations, modeling the ice
vapor interface is a particularly challenging task in terms of the quality of the force field. Ideally, the model has to adequately describe the intermolecular interactions in all three phases (solid, liquid, and vapor) and at temperatures ranging from the freezing point all the way down to ∼190 K. The presence of ions further increases the challenge, as in addition to the forces between water molecules it requires a correct description of the water
solute interactions. A large number of different potential models for water have been proposed and used over a wide range of thermodynamic conditions. For an overview of the various models and a critical comparison of their performance and predictions, we refer the reader to the extensive review by Guillot20 and the recent report on a “decathlon” test of selected water models by Vega et al.21 In the most frequently used models, water is represented as a rigid molecule with a fixed geometry, partial charges on the atoms or additional sites, and a Lennard-Jones site on oxygen to treat the intermolecular van der Waals attraction and repulsion. The most popular force fields of this kind are the SPC,20,22 SPC/E,21,23 TIP3P,22,24 and TIP4P22,24 models. These and similar models were parametrized to reproduce various properties of liquid water under ambient conditions, as they were developed primarily for simulations of biological systems. There is no guarantee, though, that a model designed for ambient conditions will be accurately transferable to other situations.20,21 Indeed, the above models often fail to reproduce the properties of gas-phase water (steam) and/or of the various solid phases (crystalline ices, amorphous solids, clathrate hydrates). Large efforts have been devoted to developing models capable of correctly describing the properties of water over a large range of phase variables. Within the framework of the simple rigid point charge models, several new potentials, including TIP5P,25,26 TIP4P-Ew,27 TIP4P/ 2005,28 TIP4P/Ice,29 and NvdE (NE6),30 have recently been designed and tested over a broad range of temperatures, pressures, and phases. In particular, the performance of these models as well as the older ones at low temperatures has been evaluated and their basic properties such as the melting temperature (see Table 1) have been calculated. Thus, over the past few years, large progress in computational studies of the structural properties of ice, the process of freezing and melting, and the properties of ice
vapor interfaces has been made.31
49 At the same time, it

became clear that there are limits to what can be achieved using rigid, nonpolarizable models.21 An obvious shortcoming of any molecular model with fixed point charges is its inability to respond, in terms of charge redistribution, to changes in the local electrostatic fields. As a consequence, the simple rigid point charge models are prone to inaccuracy in heterogeneous environments (e.g., surfaces and interfaces) or when strong polarizing forces are present in the system (e.g., due to ionic or highly polar solutes). To remedy this drawback, a number of polarizable models for water have been developed, differing in the way the polarizability of a molecule is treated, the number of sites, and other aspects (e.g., flexibility). The most recent and frequently used ones include the fluctuating charge models such as TIP4P-FQ50 or its TIP5P-based variants,51,52 the charge-on-spring (Drude oscillator) models,53,54 point dipole models (POL3,55 Dang-Chang56), the mobile charge densities in harmonic oscillators (MCDHO) model,57 and the sophisticated Thole-type model (TTM3-F).58,59 Simulations with polarizable force fields are computationally more demanding and, therefore, are much less popular. Simulation studies of ice employing polarizable water models have been particularly scarce.51,52,60
65 More importantly, unlike for the rigid, nonpolarizable models, the performance of the various polarizable models has not been extensively tested outside the narrow range of conditions for which they were originally parametrized. In particular, there are only very few polarizable water models for which the melting point has been determined. These include, to the best of our knowledge, only TIP4P-FQ (Tm = 303 K)64 and the TTM models. For the TTM3-F, the melting point calculated from classical simulations is Tm = 250 K, while Tm = 225 K was determined from quantum simulation.66 Quantum simulations were generally found to predict a lower melting temperature than their classical counterparts.66 Ultimately, it would be desirable to treat the molecular interactions in the simulated system from the first principles, without having to rely on a force field. While large progress in DFT-based molecular dynamics simulations has been made in the recent years, the computationally affordable functionals, such as Perdew
Burke
Ernzerhof (PBE) and Becke
Lee
Yang
Parr (BLYP), yield a melting temperature that is way too high (over 400 K).67 Thus, polarizable force fields currently represent an important tool with a large potential to improve our description and understanding of molecular systems. Returning back to ionic solutes, the importance of explicit treatment of polarization interactions in molecular simulations has been previously demonstrated for solvation of ions at the liquid water
vapor interface.68
72 It is reasonable to expect that polarization may be equally important also for ions at the ice
vapor interface, in particular at temperatures close to the melting point when the surface of ice is covered by a quasi-liquid layer. Simulations with nonpolarizable force fields have provided useful insights into many aspects of systems and processes involving ions and ice, such as ice nucleation in salty water,42 brine rejection from freezing salt solution,42,48,73,74 energetically favorable positions of Naþ and Cl
ions within the ice crystal lattice,75 or solvation and transfer of ions across the ice
water interface.76 However, other important questions, for example, how exactly are ionic or highly polar impurities distributed within the quasi-liquid layer on the ice surface, are likely to be adequately addressed only with polarization being explicitly taken into account. One of the popular water models that has been used extensively to characterize the solvation of ions at the surface of liquid 5974

dx.doi.org/10.1021/jp110391q |J. Phys. Chem. A 2011, 115, 5973–5982

The Journal of Physical Chemistry A water is the POL3 polarizable water potential.55 The force fields developed and used with the POL3 water model include halide anions (F
, Cl
, Br
, I
), mono- and divalent cations (Naþ, Kþ, Liþ, Mg2þ), as well as molecular ions such as nitrate (NO3
), sulfate (SO42
), hydronium (H3Oþ), and hydroxide (OH
).77
85 It would be advantageous to be able to use the same force field for investigating the behavior of ions at the suface of ice. However, as far as we know, the POL3 model has never been used for ice studies and, in particular, its melting temperature Tm is not known. We thus performed molecular dynamics simulations to determine the Tm value of the POL3 model, using the approach developed by Vega et al.,38 and to assess the suitability of this model for the simulations of ice. The results of our simulations reported in this paper show that the melting temperature of the POL3 water model is rather low (Tm = 180 ( 10 K), which severely limits its applicability to ice simulations and emphasizes the need for development of new, more reliable polarizable water models.

2. SIMULATION METHOD 2.1. Water Model. The POL3 water model55 describes a water molecule using three rigid sites located at the oxygen and hydrogen atoms, with the O
H distance of 1.0 Å and the H
O
H angle of 109.47°. The intermolecular van der Waals interaction is modeled as a Lennard-Jones potential, with a size parameter of σ = 3.2037 Å and an energy parameter of ε = 0.165 kcal/mol assigned to the oxygen site. The hydrogen sites are assigned a van der Waals term of zero. Partial charges are placed at the two hydrogens and the oxygen (qH = 0.365 e, qO =
0.73 e). In addition, each of the three sites is assigned an isotropic atomic polarizability (RH = 0.170 Å3, RO = 0.528 Å3). The polarizable POL3 model thus responds to the local environment by point dipoles μBO, μBH1, and μBH2 being induced at each of the three atom sites of a water molecule. The magnitude and direction of the induced dipoles are determined by the atomic polarizabilities RH and RO and the local electric field at the atomic positions. The electrostatic part of the intermolecular interactions then consists of two terms: the Coulomb term corresponding to the interactions between the (constant) atomic charges and the polarization term comprised of the induced dipole
charge interactions, induced dipole
induced dipole interactions, and the self-energy terms.86 While both Lennard-Jones and Coulomb interactions are pairwise additive, the polarization interaction—which is due to the response of the polarizable sites to the electric field generated by the instantaneous configuration of the atomic charges as well as by the induced dipoles on the atomic sites— inherently has a many-body character. A self-consistent routine is therefore used to evaluate the polarization part of the potential energy at every time step of the simulation. Within the POL3 model, a water molecule in a condensed phase has a permanent dipole moment (arising from the partial charges on the atoms) and an induced dipole moment (a sum of the three induced point dipoles on the polarizable atomic sites). Since the three induced point dipoles within a molecule are uncorrelated and can point in any spatial direction, the polarization response of the POL3 model is not confined to the molecular plane and does in general have an out-of-plane component. It is worth noting that the permanent molecular dipole moment of the POL3 model (2.02 D) is enhanced with respect to the experimental gas phase value (1.85 D),55 whereas the molecular polarizability (0.868 Å3) distributed between

ARTICLE

Figure 1. Initial configuration (the basal plane view) of the block of ice Ih containing 1280 water molecules. The x axis is horizontal. The ice
vacuum interfaces are located on the left- and right-hand sides of the box. The plane exposed to vacuum is the secondary prismatic (1020) plane.

oxygen and hydrogen atoms is lower than the experimental value (1.44
1.47 Å3).86 The POL3 model can therefore be regarded as a compromise between an effective pairwise potential (like SPC/ E)23 where atomic charges are enhanced to account approximately for the polarization of a water molecule in condensed phases (the average molecular dipole moment in liquid water is ∼2.6 D87) and a “pure” polarizable model where the response of a molecule to its environment is treated via polarizability rather than rescaling the charges.86 2.2. MD Simulations. To determine the melting point of the POL3 water model, we adopted the simulation protocol proposed by Vega et al.,38 employing a slab of ice Ih with two open surfaces. Another possible approach would be to determine the melting point from direct simulations of the ice
water interface31,34 or from free energy calculations.39 All of the above methods have been shown to provide melting temperatures which are in agreement with each other.38 The choice of the slab simulations has been motivated by the relative simplicity of this methodology, in particular compared to the free energy calculations. In addition, as we will show below, the melting point of the POL3 model was found to be very low. The properties of POL3 liquid water at such low temperatures are not known, and the results of the direct liquid
solid coexistence method could thus be affected by a physically unrealistic liquid phase. In all simulations, 3D periodic boundary conditions were employed and hexagonal Ih ice with 1280 water molecules in the unit cell was used. The initial configuration was constructed using the algorithm of Buch et al.44 which yields a proton disordered hexagonal ice with the overall dipole moment close to zero that fulfills the Bernal
Fowler rules.88 The dimensions of the ice sample were 44 Å  30.5 Å  29 Å. The initial configuration of ice Ih is depicted in Figure 1. To equilibrate the solid, we first performed NpT simulations of bulk ice at zero pressure38 in which the system was heated from 0 K to the desired temperature. The heating was carried out for 800
1000 ps depending on the final temperature, and the system was subsequently equilibrated at this temperature for another 200 ps. In comparison with the rigid nonpolarizable water models,38 the POL3 ice structure turned out to be rather sensitive during 5975

dx.doi.org/10.1021/jp110391q |J. Phys. Chem. A 2011, 115, 5973–5982

The Journal of Physical Chemistry A annealing. Therefore, a very short time step of 0.1 fs was used during the heating and equilibration phase, in addition to a slow increase of temperature. Large disordering or even melting of the ice structure often occurred when a longer time step and/or faster heating were used. It should be noted, however, that it is not uncommon for a model ice structure to exhibit instability during annealing. While heating of the POL3 ice was a particularly delicate procedure, such behavior is not an exception specific to this particular force field. Heating of ice is a process which has to be carried out with great care for the nonpolarizable models, too, to avoid large angular reorientations of water molecules that can trigger collective motions leading to disordering of the ice structure during the heating phase.89 The temperature was regulated using the Berendsen thermostat90 with a relaxation time of 0.1 ps. The constant pressure was maintained by a barostat employing the weak-coupling scheme analogous to the temperature control algorithm with a relaxation time of 2 ps. The pressure of the barostat was set to zero. All three sides of the simulation box were allowed to fluctuate independently in order not to a priori impose the dimensions of the unit cell. Rather, the equilibrium size of the simulation box is defined by the system, i.e., by the interaction potential and the thermodynamic conditions. To evaluate the long-range electrostatic interactions, particle mesh Ewald summation91,92 was used. For Lennard-Jones interactions and the real space part of the electrostatic interaction, a cutoff distance of 11 Å was applied. Water bond lengths and angles were constrained using the SHAKE algorithm.93 Upon equilibration of the bulk ice (at zero pressure) to the temperature of interest, the ice slab was prepared by extending one dimension of the box (in our case the x axis) from its original value of ∼44 to 100 Å. In this way, an approximately 44 Å thick ice slab with two ice
vacuum interfaces was created (the ice slab thickness of ∼40 Å has been shown to yield a stable block of ice that melts at the melting temperature38). As a result, the size of the box in the slab simulations was 100 Å  ∼30.5 Å  ∼29 Å, with the x axis perpendicular to the ice
vacuum interface. The exact lateral (y and z) dimensions of the simulation box were determined by the preceding NpT run at each particular temperature. The periodic images of the slab in the x direction were separated by an ∼56 Å thick layer of vacuum (vapor phase). The plane which was exposed to the vacuum was a secondary prismatic (1020) plane. For more details regarding the main crystallographic planes of ice (basal, primary prismatic, and secondary prismatic), see for instance the water Web site of Chaplin.94 The selection of the plane to form the ice
vacuum interface is arbitrary; despite the fact that the properties of the ice surface and, in particular, the melting mechanism may depend on the crystal plane facing the vacuum, the melting temperature should be independent of this choice. However, the selection of the secondary prismatic plane has its advantages, as it has been shown for the ice
water interfaces that the secondary prismatic plane exhibits the fastest dynamics.49,95 In addition, choice of this face as the free surface is convenient for the visualization of the degree of order/disorder during the ice melting process. Once the simulation box was extended and the two ice
vacuum interfaces were created, we proceeded with a set of relatively long NVT simulations (up to 30 ns) of the ice slab at several different temperatures to determine the lowest temperature at which the ice slab melts, Tþ, and the highest temperature at which the slab remains stable and does not melt, T
. The melting point Tm of the model is then estimated as Tm = 1/2(Tþ

ARTICLE

Figure 2. Total energy per one molecule as a function of time in MD simulations of an ice Ih slab at various temperatures. The time t = 0 corresponds to the moment when the ice slab was created, i.e., when the simulation box containing pre-equilibrated bulk ice was extended along the x axis, exposing thus the secondary prismatic plane to the vapor phase.

þ T
).38 In the slab simulations (NVT), the dimensions of the box were fixed and no pressure control was used. A time step of 1 fs was employed. The motion of the center of mass of the system was removed every picosecond. All the other conditions and control parameters of the slab simulations were identical to the ones used in the bulk ice (NpT) simulations. The data in the slab simulations were recorded every 10 ps for analysis from the moment the configuration was changed from bulk to slab geometry. In this way, the changes in energy related to the formation of the ice
vapor interfacial region could be captured. In addition, a control simulation of bulk liquid water at T = 180 K was performed in order to obtain an independent estimate of the enthalpy of melting of the POL3 model. The number of water molecules (1280) in the liquid water simulation was the same as that used in the ice sample. A box of liquid water equilibrated at p = 1 bar and T = 300 K was used as the initial condition. The system was then simulated in an NpT run at zero pressure and T = 180 K. The total energy of the liquid was compared with the bulk ice simulated under the same conditions to determine the energy difference between the ice and liquid water at the melting temperature. All simulations were performed using the AMBER 8.0 program package.96 Images of the simulated systems were created using the VMD (Visual Molecular Dynamics) program.97

3. RESULTS In Figure 2, the results for total energy per one molecule of the ice slab throughout the MD runs at different temperatures are presented. In addition, density profiles averaged over the last 4 ns of the trajectories and final snapshots of the simulated slabs for a selected subset of temperatures are shown in Figures 3 and 4. At the highest simulated temperature (210 K), there is an increase of the total energy at the beginning of the simulation (first 3
4 ns) and then the energy levels off and remains approximately constant, subject to fluctuations. Similar behavior is observed at 200 K, albeit the initial increase in energy is less steep and is extended over a somewhat longer time (∼6 ns). The visualization of the trajectories as well as the corresponding density profiles indicate that the ice 5976

dx.doi.org/10.1021/jp110391q |J. Phys. Chem. A 2011, 115, 5973–5982

The Journal of Physical Chemistry A

ARTICLE

Figure 3. Density profiles of water oxygen atoms of the slab system for selected temperatures averaged over the last 4 ns of each simulation. The position, r, is measured across the slab, in the direction normal to the ice
vapor interface.

Figure 4. Instantaneous configurations of the system at the end of the simulation at selected temperatures (total length of each simulation given in parentheses).

slab melts (see Figures 3a and 4a). However, even though the system clearly is above its melting point and only a few nanoseconds are needed for the ice slab to melt (∼4 ns at 210 K and ∼6 ns at 200 K), a much longer time is required afterward to obtain fully equilibrated liquid. The diffusion of water molecules in the liquid

at ∼200 K is very slow, and as can be seen in the density profile (Figure 3a), traces of the original crystalline structure persist in the melted system over a relatively long time scale (>10 ns). With decreasing the temperature to 195 K and below, the melting process is further slowed down and the total energy 5977

dx.doi.org/10.1021/jp110391q |J. Phys. Chem. A 2011, 115, 5973–5982

The Journal of Physical Chemistry A remains almost constant during the entire simulation (up to ∼30 ns); only thermal fluctuations with a minimal (if any) overall increase in energy can be seen. Nevertheless, from the visual inspection of the system, it can be seen on the simulated time scale that down to 185 K the ice slab melts. The final snapshot of the system at 185 K at the end of the 27 ns simulation is shown in Figure 4b. Melting of the ice slab is confirmed by the corresponding density profile shown in Figure 3b. However, both the density profile and the snapshot show that a large degree of crystalline order persists in the melting system over the time scale of tens of nanoseconds. The dynamics of melting is extremely slow at these temperatures, so that it may take an immensely long time for the system to melt completely. Regardless of the incomplete melting, though, the ice slab is clearly not stable. Let us now turn our attention to low temperatures. At 160
180 K, rapid formation of a quasi-liquid layer on both open surfaces of the ice slab is observed, while the middle part of the slab remains solid throughout the simulation (see the density profiles in Figure 3c,d and the snapshots in Figure 4c,d). It can be observed that the thickness of the quasi-liquid layer is spatially rather nonuniform, though on average it decreases with decreasing temperature as expected. In addition, as can be seen, the bulk ice part of the slab exhibits a certain degree of disorder even at the lowest simulated temperature (160 K), and this disorder is increasing with increasing temperature. On the basis of the above results, we may conclude that the melting point of the POL3 model is located below 185 K. However, the extremely slow dynamics of the melting process around this temperature makes the assignment of Tþ (lowest temperature at which the ice slab melts) and T
(highest temperature at which the slab does not melt) very difficult. Tþ and T
cannot be easily determined from the plot of the total energy as a function of time at different temperatures (Figure 2) as in a broad range of temperatures around the melting point the energy remains approximately constant and there is no clear difference between the time evolution of the energy in the case when the slab melts or when it remains frozen. The most straightforward criterion to establish whether the ice slab is melting or not seems to be the visualization of the system combined with the monitoring of the density profile of the slab. Unfortunately, the results of this approach are dependent on the time scale of the MD simulations. The present set of simulations indicates that the melting point of the POL3 model may be located between 180 K, where a portion of crystalline ice is clearly visible in the interior of the slab, and 185 K, where the ice slab clearly melts. However, despite the fact that this assignment of T
and Tþ is based on rather long simulation times (up to ∼30 ns), we cannot rule out the possibility that over a substantially longer time scale the ice slab would melt at 180 K (or even at a lower temperature). Therefore, we prefer to make a somewhat more conservative estimate of the melting temperature and we assign the melting point of the POL3 water to be Tm = 180 ( 10 K. Further refinement of the melting point value and of the error associated with Tm (in other words, the difference between Tþ and T
) can in principle be achieved by performing simulations at additional, more densely spaced temperatures and/or by carrying out longer MD runs (of the order of hundreds rather than tens of nanoseconds). However, taking into account that the diffusion of the POL3 water close to the melting point is extremely slow and, in addition, the simulations with explicit treatment of polarization are computationally costly, the above

ARTICLE

Table 2. Average Values of the Dipole Moment and the Components of the Quadrupole Moment of a Water Molecule in Liquid Water and in Ice (at T = 180 K) Modeled by the POL3 Water Modela μ (D) Qxx (DÅ) Qyy (DÅ) Qzz (DÅ) QT (DÅ) μ/QT (Å
1) permanent 2.026

1.89


1.62

induced

0.727

0.14

0

total

2.736

2.03


1.62

induced

0.714

0.16

0

total

2.726

2.05


1.62


0.26

1.76

1.154

1.83

1.169

1.84

1.168

Liquid 0.11
0.15

Ice 0.11
0.16

QT is defined as QT = (Qxx
Qyy)/2. The origin is located in the center of mass of the water molecule, the z axis coincides with the H
O
H angle bisector, and the x axis is parallel to the vector connecting the two hydrogen atoms. a

result represents what could be accomplished with our computer resources within reasonable time. In their paper, Abascal and Vega36 conclude that the melting point is strongly dependent on the quadrupole moment of the water models, while it seems almost uncorrelated with their dipole moment. To investigate this aspect for the POL3 model, we evaluated the average molecular dipole and quadrupole moments for liquid water and ice at 180 K. Following refs 21 and 61, we chose the origin for the quadrupole moment calculation to be located in the center of mass of the water molecule. The z axis coincides with the H
O
H angle bisector, the x axis is parallel to the vector connecting the two hydrogen atoms, and the y axis is normal to the molecular plane. The results are shown in Table 2. The permanent dipole and quadrupole moments of the POL3 water molecule are given by the rigid molecular geometry (the O
H distance of 1.0 Å and the H
O
H angle of 109.47°) and the fixed partial charges located at the atom sites (qH = 0.365 e, qO =
0.73 e). For the above choice of the molecular axes, the off-diagonal elements of the permanent quadrupole tensor vanish.36 The induced dipole of a water molecule is given by a sum of the three induced point dipoles on the polarizable atomic sites. The total molecular dipole moment is, in turn, obtained as a vector sum of the permanent and induced dipoles. As seen in Table 2, for the POL3 model, the average induced water dipole amounts to about one-third of the permanent dipole both in liquid water and in ice at 180 K (i.e., T ∼ Tm). The average total molecular dipole is about 2.7 D in both cases, which is somewhat higher than the value of 2.6 D reported for POL3 water at ambient temperature.86,98 Note that, while the permanent water dipole is confined to the molecular plane, the induced dipole can point in any direction. Therefore, the magnitude of the total dipole is not in general equal to the sum of the magnitudes of the permanent and induced dipole moments. The out-of-plane component of the induced dipoles was, however, found to be rather small in our simulations, amounting on average to ∼0.1 D. In addition to the permanent quadrupole due to the atomic charges, the point dipoles μBO, μBH1, and μBH2 at the polarizable atomic sites give rise to an induced quadrupole moment given by ƒ! X r þ μ ƒ! X r Q5 ¼ ƒ! μO X B rOþμ H1 H2 BH1 BH2

ð1Þ

whereBr O,Br H1, andBr H2 are the position vectors of the oxygen and 5978

dx.doi.org/10.1021/jp110391q |J. Phys. Chem. A 2011, 115, 5973–5982

The Journal of Physical Chemistry A the two hydrogen atoms with respect to the center of mass of the water molecule and X denotes the tensor product of the two respective vectors. In our coordinate system, all three position vectors lie within the molecular plane, hence Qyy = 0. The induced quadrupole tensor in principle has nonzero off-diagonal components; however, analysis of our simulation data reveals that the average values of the off-diagonal elements are at least 4 orders of magnitude smaller than the diagonal elements Qxx and Qzz. At the same time, these two nonzero diagonal elements of the induced quadrupole tensor are substantially smaller than the corresponding elements of the permanent quadrupole (see Table 2). Thus, the induced quadrupole represents only a minor contribution to the total quadrupole of POL3 water, and the properties of the model are determined primarily by the permanent quadrupole due to the fixed charge distribution within the molecule. This allows us to use, albeit approximately, the quantity QT = (Qxx
Qyy)/2 as a significant measure of the magnitude of the total molecular quadrupole moment for the purpose of comparison of POL3 with the nonpolarizable models, as discussed below.

4. DISCUSSION It is interesting to compare the present results for the POL3 ice slab to those obtained by Vega et al. for nonpolarizable water models (SPC/E, TIP4P, TIP4P/Ew, TIP4P/2005, TIP4P/ Ice).38 Contrary to the melting behavior of the nonpolarizable models, the total energy curves of the POL3 model in the simulations above Tm do not exhibit a final steep increase in energy followed by a plateau corresponding to completed melting and formation of liquid in equilibrium with its vapor. We note that the melting curves of the POL3 model most closely resemble those of the SPC/E water model, that is, the three-site model upon which the POL3 model was constructed. Out of the models studied in ref 38, the SPC/E is the only three-site model and has by far the lowest melting temperature of ∼215 K. (The other four are four-site models of the TIPnP family with melting points between 230 and 271 K.) In comparison to the four-site models, the time evolution of the total energy of the SPC/E model above its melting point (cf. Figure 6 of ref 38) exhibits a long (up to ∼15 ns) period of slow increase of energy. Apparently, in the case of the three-site SPC/E model with its very low melting point, longer times are needed to achieve complete melting of the ice sample. Melting of the POL3 ice slab in the present simulations proceeds even slower. The POL3 model differs from the nonpolarizable water models also as far as the enthalpy of melting is concerned. The difference in total energy per molecule found in the present study between the beginning of the simulation (frozen slab) and the end of the simulation (melted slab) is about 0.7 kJ/mol at 210 K and about 0.5 kJ/mol at 200 K. We performed a control NpT run of liquid water at p = 0 bar and T = 180 K and compared the total energy of this system to bulk ice simulated under the same conditions. The total energy difference of 0.75 kJ/mol was found, corresponding well to the above values observed between the frozen and melted ice slab. Apparently, the melting enthalpy of the POL3 model is rather low compared to the experimental value of 6.03 kJ/mol as well as to the melting enthalpies of about 2
4 kJ/mol reported for the nonpolarizable water models.38 For this reason, there is only a very small increase of energy upon melting of the ice slab (Figure 2).

ARTICLE

Another difference between the POL3 ice and ice modeled by nonpolarizable water force fields concerns the degree of disorder present in the bulk ice. In addition to the well-known L- and D-defects introduced by Bjerrum, vacancies and similar nonlocal defects, there are also defects caused by local rearrangement of hydrogen bonds in which six-membered rings around a hydrogen bond are replaced by 5- and 7-rings (the so-called “5 þ 7 defect”).99 While recent quantum chemical calculations100 suggest that the defect is real and its probability is around 4  10
8 to 2  10
9 (per water molecule), the probabilities of the defect predicted by empirical force fields vary by several orders of magnitude depending on the water model used to construct the ice. It has thus been proposed that the hydrogen-bond defects in the structure of ice may be used as a test of quality of empirical force fields employed for simulations of water ice Ih. From visualization of the POL3 ice slab at different temperatures in the present simulations (see Figure 4c,d), it is clear that the degree of hydrogen-bond disorder occurring in the interior of the block of POL3 ice is substantially larger (at the corresponding degree of undercooling) than in ice modeled by nonpolarizable water models, including the TIP4P/2005 and TIP4P/Ice foursite models38 (which predict the probability of the 5 þ 7 defect about 10
5)100 as well as the three-site ones such as SPC/E.73 It should be noted that the increased disorder present in the POL3 ice is not specific to the slab geometry (or an artifact thereof). We have visualized and analyzed all runs from their very beginning, starting with the heating of the bulk ice from 0 K to the desired temperature and subsequent equilibration at the target temperature, before a slab of ice was created by extending the simulation box along the dimension normal to the ice
vapor interface. The higher level of hydrogen-bond disorder relative to the nonpolarizable water models was found to be the case also in the bulk ice and was observed already at rather low temperatures during the heating of the ice. Consequently, as discussed above in the Methods section, heating of the POL3 ice required extra care and very slow heating rate. The above observations highlight the importance of adequate slab thickness in MD simulations of ice
vapor interface. It has been concluded from previous simulations38 employing nonpolarizable water models that the thickness of the block of ice has to be sufficiently large, not smaller than about 3
4 nm, in order to achieve a slab with two stable ice
vapor interfaces which melts at the melting temperature. For a stable ice slab, the melting process starts at the surface and then, for T > Tm, propagates to the interior of the slab, whereas for T < Tm the melting is limited to the surface layer (QLL) while the interior of the slab remains crystalline. If the layer of ice is not thick enough and the two ice
vapor interfaces are not well separated, the surface premelting may propagate from both interfaces to the interior of the ice slab and result in complete melting of the block of ice even at T < Tm. The present study indicates that for the POL3 model larger slab thickness (g4 nm) is advisible in order to avoid such instability. It should also be mentioned that the planar geometry of the POL3 model may influence the rotational vs translational disordering around the melting point. Indeed, as can be seen in Figure 4, the rotational disordering appears to be more pronounced at melting, whereas the translational order persists up to 200 K (Figure 3a). It has been shown that, near the melting point, the orientational disorder and the corresponding rotational mobility are larger than the translational ordering and mobility for planar water models,101,102 in contrast with the nonplanar 5979

dx.doi.org/10.1021/jp110391q |J. Phys. Chem. A 2011, 115, 5973–5982

The Journal of Physical Chemistry A TIP5P water potential.102 The persistence of the translational order at T > Tm may thus highlight unphysical description of the liquid phase at these low temperatures by the POL3 potential. As regards the melting point of the POL3 model, a rather low value of Tm was to be expected on the basis of the fact that threesite models in general yield very low melting temperature, while models with more than three sites representing the water molecule are closer to experiment (see Table 1). As has been shown by Abascal and Vega,36 the failure of the three-site models to yield correct melting temperature is a consequence of their low quadrupole moments which, in turn, is a consequence of the negatively charged site being located at the oxygen atom. In particular, the geometry of the water molecule (i.e., the bond length and angle) and the distribution of the permanent charge sites in the POL3 model is the same as in SPC and SPC/E, which both exhibit melting around 200 K. The polarization correction,103 which transformed SPC (qH = 0.41e) into SPC/E (qH = 0.423e) by taking into account the averaged increase of the molecular dipole moment of water (via enlarging the atomic partial charges) in the liquid phase and which resulted in improved liquid properties, yielded also a somewhat better melting temperature (SPC, Tm = 190 K; SPC/E, Tm = 215 K). In POL3, the charges on atomic sites are decreased again (qH = 0.365e), yielding a lower permanent dipole moment (2.02 D) relative to that of SPC/E (2.35 D) and SPC (2.27 D). An additional dipole moment due to the induction forces in the liquid phase is calculated explicitly by introducing polarizable centers on the atomic sites, resulting in an averaged total dipole moment of a POL3 water molecule86,98 experimental value of 2.6 D.87,104 While this is certainly a refinement with respect to the nonpolarizable models as far as the treatment of the many-body polarization forces are concerned, one might expect Tm of the POL3 water to be comparable if not better to Tm of the SPC/E model. However, the performance of the POL3 force field in terms of the melting point is actually worse than that of the original SPC model which does not even account for the condensed-phase polarization in an effective way via the enhanced charges. This can be rationalized by taking into account also the quadrupolar interactions, the importance of which for the properties of the water solid phases, including the melting point of ice Ih, has recently been demonstrated.21 Our results for the POL3 model show that, while the induced point dipoles at the polarizable atomic centers contribute substantially to the total molecular dipole (the average size of the induced dipole on a water molecule is about one-third of the permanent dipole moment), the contribution of the same induced point dipoles to the total quadrupole moment is rather small. For POL3 liquid water at 300 K, even a small net reduction of the quadrupole moment due to the induced dipoles was reported in the original work of Caldwell and Kollman.55 Thus, the properties of the model are determined primarily by the permanent quadrupole, which is by itself very low due to both the three-site geometry and rather low values of atomic charges. Using the quantity QT = (Qxx
Qyy)/2 as a measure of the magnitude of the molecular quadrupole moment21 reveals QT = 1.76 DÅ for the permanent quadrupole and QT = 1.83 (1.84) DÅ for the total quadrupole in water (ice) in our simulations at T ∼ Tm. Comparison of the above value of QT for the total molecular quadrupole against the nonpolarizable models (see Table 2 in ref 21) places the POL3 model below both SPC/E (QT = 2.035 DÅ) and SPC (QT = 1.969 DÅ) while still above the TIP3P model (QT = 1.721 DÅ). The same holds for the melting points of these models. In

ARTICLE

addition, as can be seen in Table 2, the ratio of the dipole to quadrupole moment, μ/QT, which seems to play a crucial role in the quality of the phase diagram predicted by the different water models,21 deviates substantially from 1 for the POL3 model. A poor description of the phase diagram can, therefore, be expected as a result of disbalance between dipolar and quadrupolar forces. Indeed, the fact that the melting point of ice Ih is so low when using the POL3 model strongly suggests that ice II or another ice polymorph is the stable phase at melting, as has been shown, for instance, for the TIP3P model.21 As for the TIP3P model, one may anticipate that the melting temperature of the stable ice phase for the POL3 model is above that of the Ih phase determined here. While a detailed characterization of the phase diagram of POL3 ice is beyond the scope of the present study, the results presented here clearly indicate that this model is rather unsuitable for ice simulations.

5. CONCLUSIONS Molecular dynamics simulations of a slab of ice Ih modeled using the polarizable POL3 force field at a set of different temperatures were performed following the approach of Vega et al.38 to determine the melting temperature of the POL3 model and to assess the suitability of this model for the simulations of ice. The simulations yielded the melting point of the POL3 model to be Tm = 180 ( 10 K. While it has been shown that quadrupolar interactions play a critical role in the melting temperatures of simple (nonpolarizable) water models and acceptable predictions for Tm can only be obtained when the negatively charged site is shifted away from the oxygen atom toward the hydrogens (which is not the case in the three-site models including POL3), the melting temperature of the POL3 model turned out to be significantly lower than Tm of the SPC/E model and even ∼10° below Tm of the original SPC force field which, unlike SPC/E, does not include the effective polarization correction. This is rationalized in terms of rather different contributions of the induced point dipoles at the three atomic polarizable sites of the POL3 model to the total molecular dipole as opposed to the total molecular quadrupole. While the induced point dipoles contribute substantially to the overall polarization of the water molecule and the total molecular dipole moment increases significantly in a condensed phase, the induced quadrupole moment is very small and, thus, the total molecular quadrupole in a condensed phase remains small as a result of the three-site molecular geometry combined with rather low atomic charges. This disbalance between the dipolar and quadrupolar interactions is responsible for the low melting point of ice Ih and indicates that the POL3 model is likely to yield a poor description of the phase diagram. In particular, the fact that the melting point is so low strongly suggests that ice Ih does not constitute a stable phase at temperatures around Tm. Thus, the results clearly show that the POL3 water force field is rather a poor model for studying ice and ice
liquid or ice
vapor interfaces. A straightforward conclusion would be to avoid the use of the POL3 model for such studies. The situation is, however, not that simple. The number of reliable polarizable force fields that are presently available for ice simulations is very limited. While the pairwise additive, nonpolarizable models often provide a sufficient, though not ideal description of the ice environment, such approximation may be too severe in systems containing ions or highly polar molecules in which induction forces are significant. Moreover, different 5980

dx.doi.org/10.1021/jp110391q |J. Phys. Chem. A 2011, 115, 5973–5982

The Journal of Physical Chemistry A polarizable force fields treat the many-body polarization interaction in different ways. Therefore, the force field parametrizations, developed for a solute in conjunction with a certain water model, often cannot be used with another water model, and reparametrization or development of an entirely new force field for the solute is necessary. Over the past years, force fields for a large number of ions and other solutes have been parametrized and successfully used in simulations with the POL3 water model. Thus, in situations when polarizability is crucial (for example, when the detailed solvation structure of large “soft” ions such as Br
or I
or ions with nonspherical molecular structure such as NO3
plays an important role) and use of a polarizable force field is necessary, POL3 may still at present be the model of choice and, if employed with caution, can provide useful results despite its severe drawbacks. Further work is necessary to evaluate the low-temperature properties of the POL3 model in more detail, in particular as regards the ice
liquid coexistence and the quasiliquid layer at the ice
vapor interface, because liquid water at T ∼ Tm = 180 K is deep below the ambient conditions for which the model was developed, and this may severely affect both its structure and dynamics. In summary, this study shows severe deficiencies of the POL3 water model when applied to the simulations of ice and ice
vapor interface systems. At the same time, it highlights the urgent need for a new polarizable water force field, reliable over larger areas of the phase diagram. Since the three-site models are clearly not able to yield correctly even such fundamental quantity as the melting temperature, inclusion of polarizability within the water models with four or more sites should be explored as a possible way to further improvement.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT We would like to thank Carlos Vega (Universidad Complutense, Madrid), Jirí Kolafa (Institute of Chemical Technology, Prague), and Pavel Jungwirth (Institute of Organic Chemistry and Biochemistry, Prague) for helpful discussions. We also thank two anonymous reviewers for their insightful comments. This work has been supported by the Czech Science Foundation (grant P208/10/1724), the Ministry of Education of the Czech Republic (grants MEB020919 and LC512), and the Academy of Sciences of the Czech Republic (project AV0Z40550506). Support from the bilateral French/Czech BARRANDE program under Project No. 14248QK (French side)/MEB020715 (Czech side) and by the French Ministry of Research through the LEFE/ CHAT program is also gratefully acknowledged. ’ REFERENCES (1) Barker, J. A.; Watts, R. O. Chem. Phys. Lett. 1969, 3, 144–145. (2) Rahman, A.; Stillinger, F. H. J. Chem. Phys. 1971, 55, 3336–3359. (3) Heymsfeld, A. J.; Sabin, R. M. J. Atmos. Sci. 1989, 46, 2252. (4) Hoyle, C. R.; Luo, B. P.; Peter, T. J. Atmos. Sci. 2005, 7, 2568. (5) Sassen, K.; Liou, K. N.; Kinnes, S.; Griffin, M. Science 1985, 227, 410. (6) Seinfeld, J. H.; Pandis, S. N. Atmos. Chem. Phys.; Wiley: New York, 1998. (7) K€archer, B.; Voigt, C. Geophys. Res. Lett. 2006, 33, L08806.

ARTICLE

(8) Kr€amer, M.; Schiller, C.; Ziereis, H.; Overlez, J.; Bunz, H. Tellus, Ser. B 2006, 58B, 141. (9) Stuart, A. L.; Jacobson, M. Z. J. Geophys. Res. 2004, 108, 425. (10) Collignon, B.; Picaud, S. Chem. Phys. Lett. 2004, 393, 457–463. (11) Peybernes, N.; Le Calve, S.; Mirabel, P.; Picaud, S.; Hoang, P. N. M. J. Phys. Chem. B 2004, 108, 17425–17432. (12) Picaud, S.; Hoang, P. N. M. J. Chem. Phys. 2000, 112, 9898–9908. (13) Picaud, S.; Hoang, P. N. M.; Peybernes, N.; Le Calve, S.; Mirabel, P. J. Chem. Phys. 2005, 122. (14) Petitjean, M.; Hantal, G.; Chauvin, C.; Mirabel, P.; Le Calv
e, S.; Hoang, P. N. M.; Picaud, S.; Jedlovszky, P. Langmuir 2010, 26, 9596. (15) Hantal, G.; Jedlovszky, P.; Hoang, P. N. M.; Picaud, S. J. Phys. Chem. C 2007, 111, 14170. (16) Hantal, G.; Jedlovszky, P.; Hoang, P. N. M.; Picaud, S. Phys. Chem. Chem. Phys. 2008, 10, 6369. (17) Jedlovszky, P.; Hantal, G.; Neurohr, K.; Picaud, S.; Hoang, P. N. M.; Von Hessberg, P.; Crowley, J. N. J. Phys. Chem. C 2008, 112, 8976. (18) Jedlovszky, P.; Partay, L.; Hoang, P. N. M.; Picaud, S.; von Hessberg, P.; Crowley, J. N. J. Am. Chem. Soc. 2006, 128, 15300. (19) P
artay, L. B.; Jedlovszky, P.; Hoang, P. N. M.; Picaud, S.; Mezei, M. J. Phys. Chem. C 2007, 111, 9407. (20) Guillot, B. J. Mol. Liq. 2002, 101, 219–260. (21) Vega, C.; Abascal, J. L. F.; Conde, M. M.; Aragones, J. L. Faraday Discuss. 2009, 141, 251–276. (22) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Intermolecular forces; Reidel: Dordrecht, The Netherlands, 1981. (23) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (24) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (25) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (26) Rick, S. W. J. Chem. Phys. 2004, 120, 6085–6093. (27) Horn, H. W.; Swope, W. C.; Pitera, J. W.; Madura, J. D.; Dick, T. J.; Hura, G. L.; Head-Gordon, T. J. Chem. Phys. 2004, 120, 9665. (28) Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2005, 123, 234505. (29) Abascal, J. L. F.; Sanz, E.; Garcia Fern
andez, R.; Vega, C. J. Chem. Phys. 2005, 122, 234511. (30) Nada, H.; van der Eerden, J. P. J. M. J. Chem. Phys. 2003, 118, 7401. (31) Wang, J.; Yoo, S.; Bai, J.; Morris, J. R.; Zeng, X. C. J. Chem. Phys. 2005, 123, 036101. (32) Vrbka, L.; Jungwirth, P. J. Phys. Chem. B 2006, 110, 18126. (33) Bolton, K.; Petterson, J. B. C. J. Phys. Chem. B 2000, 104, 1590–1595. (34) Garcia Fernandez, R.; Abascal, J. L. F.; Vega, C. J. Chem. Phys. 2006, 124, 144506. (35) Abascal, J. L. F.; Garcia Fernandez, R.; Vega, C.; Carignano, M. A. J. Chem. Phys. 2006, 125, 166101. (36) Abascal, J. L. F.; Vega, C. Phys. Chem. Chem. Phys. 2007, 9, 2775–2778. (37) Conde, M. M.; Vega, C.; Patrykiejew, A. J. Chem. Phys. 2008, 129, 014702. (38) Vega, C.; Martin-Conde, M.; Patrykiejew, A. Mol. Phys. 2006, 104, 3583–3592. (39) Vega, C.; Sanz, E.; Abascal, J. L. F. J. Chem. Phys. 2005, 122, 114507. (40) Bishop, C. L.; Pan, D.; Liu, L. M.; Tribello, G. A.; Michaelides, A.; Wang, E. G.; Slater, B. Faraday Discuss. 2009, 141, 277–292. (41) Neshyba, S.; Nugent, E.; Roeselov
a, M.; Jungwirth, P. J. Phys. Chem. C 2009, 113, 4597–4604. (42) Bauerecker, S.; Ulbig, P.; Buch, V.; Vrbka, L.; Jungwirth, P. J. Phys. Chem. C 2008, 112, 7631–7636. (43) Buch, V.; Sadlej, J.; Aytemiz-Uras, N.; Devlin, J. P. J. Phys. Chem. A 2002, 106, 9374–9389. 5981

dx.doi.org/10.1021/jp110391q |J. Phys. Chem. A 2011, 115, 5973–5982

The Journal of Physical Chemistry A (44) Buch, V.; Sandler, P.; Sadlej, J. J. Phys. Chem. B 1998, 102, 8641–8653. (45) Cwiklik, L.; Devlin, J. P.; Buch, V. J. Phys. Chem. A 2009, 113, 7482–7490. (46) Devlin, J. P.; Buch, V. J. Chem. Phys. 2007, 127, 091101. (47) Groenzin, H.; Li, I.; Buch, V.; Shultz, M. J. J. Chem. Phys. 2007, 127, 214502. (48) Carignano, M. A.; Baskaran, E.; Shepshon, P. B.; Szleifer, I. Ann. Glaciol. 2006, 44, 113. (49) Carignano, M. A.; Shepshon, P. B.; Szleifer, I. Mol. Phys. 2005, 103, 2957–2967. (50) Rick, S. W.; Stuart, S. J.; Berne, B. J. J. Chem. Phys. 1994, 101, 6141. (51) Baranyai, A. J. Mol. Liq. 2009, 148, 88–93. (52) Baranyai, A.; Bartok, A. J. Chem. Phys. 2007, 126, 184508. (53) Yu, H.; Hansson, T.; van Gunsteren, W. F. J. Chem. Phys. 2003, 118, 221. (54) Lamoureux, G.; MacKerell, J. A. D.; Roux, B. J. Chem. Phys. 2003, 119, 5185. (55) Caldwell, J. W.; Kollman, P. A. J. Phys. Chem. 1995, 99, 6208–6219. (56) Dang, L. X.; Chang, T.-M. J. Chem. Phys. 1997, 106, 6772. (57) Saint-Martin, H.; Hern
andez-Cobos, J.; Bernal-Uruchurtu, M. I.; Ortega-Blake, I.; Berendsen, H. J. C. J. Chem. Phys. 2000, 113, 10899. (58) Fanourgakis, G. S.; Xantheas, S. S. J. Chem. Phys. 2008, 128, 074506. (59) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 5115. (60) Saint-Martin, H.; Hess, B.; Berendsen, H. J. C. J. Chem. Phys. 2004, 120, 11133–11143. (61) Rick, S. W. J. Chem. Phys. 2001, 114, 2276–2283. (62) Rick, S. W. J. Chem. Phys. 2005, 122, 094504. (63) Rick, S. W.; Haymet, A. D. J. J. Chem. Phys. 2003, 118, 9291. (64) Nicholson, B. F.; Clancy, P.; Rick, S. W. J. Cryst. Growth 2006, 293, 78–85. (65) Fanourgakis, G. S.; Xantheas, S. S. J. Chem. Phys. 2006, 124, 174504–174504. (66) Paesani, F.; Yoo, S.; Bakker, H. J.; Xantheas, S. S. J. Phys. Chem. Lett. 2010, 1, 2316–2321. (67) Yoo, S.; Zeng, X. C.; Xantheas, S. S. J. Chem. Phys. 2009, 130, 221102. (68) Dang, L. X. J. Phys. Chem. B 2002, 106, 10388. (69) Dang, L. X. J. Chem. Phys. 2003, 119, 6351. (70) Dang, L. X.; Chang, T. S.; Roeselova, M.; Garrett, B. C.; Tobias, D. J. J. Chem. Phys. 2006, 124, 066101. (71) Jungwirth, P.; Tobias, D. J. Chem. Rev. 2006, 106, 1259–1281. (72) Mucha, M.; Frigato, T.; Levering, L. M.; Allen, H. C.; Tobias, D. J.; Dang, L. X.; Jungwirth, P. J. Phys. Chem. B 2005, 109, 7617–7623. (73) Vrbka, L.; Jungwirth, P. Phys. Rev. Lett. 2005, 95, 148501. (74) Carignano, M. A.; Shepshon, P. B.; Szleifer, I. Chem. Phys. Lett. 2007, 436, 99–103. (75) Smith, E. J.; Haymet, A. D. J. Mol. Simul. 2004, 30, 827–830. (76) Smith, E. J.; Bryk, T.; Haymet, A. D. J. J. Chem. Phys. 2005, 123, 034706. (77) Dang, L. X. J. Chem. Phys. 1992, 96, 6970. (78) Dang, L. X.; Garrett, B. C. J. Chem. Phys. 1993, 99, 2972– 2977. (79) Cannon, W. R.; Pettitt, B. M.; McCammon, J. A. J. Phys. Chem. 1994, 98, 6225. (80) Salvador, P.; Curtis, J. E.; Tobias, D. J.; Jungwirth, P. Phys. Chem. Chem. Phys. 2003, 5, 3752–3757. (81) Wahab, A.; Mahiuddin, S.; Hefter, G.; Kunz, W.; Minofar, B.; Jungwirth, P. J. Phys. Chem. B 2005, 109, 24108–24120. (82) Perera, L.; Berkowitz, M. L. J. Chem. Phys. 1994, 100. (83) Callahan, K.; Casillas-Ituarte, N. N.; Roeselova, M.; Allen, H. C.; Tobias, D. J. J. Phys. Chem. A 2010, 114, 5141–5148. (84) Vacha, R.; Megyes, T.; Bako, I.; Pusztai, I.; Jungwirth, P. J. Phys. Chem. A 2009, 113, 4022.

ARTICLE

(85) Buch, V.; Milet, A.; Vacha, R.; Jungwirth, P.; Devlin, J. P. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 7342. (86) Kolafa, J. Collect. Czech. Chem. Commun. 2008, 73, 507–517. (87) Coulson, C.; Eisenberg, D. Proc. R. Soc. London, Ser. A 1966, 291, 445. (88) Bernal, J. D.; Fowler, R. H. J. Chem. Phys. 1933, 1, 515. (89) Girardet, C.; Toubin, C. Surf. Sci. Rep. 2001, 44, 163–238. (90) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (91) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (92) Darden, T.; York, D.; Pedersen, L. G. J. Chem. Phys. 1993, 98, 10089. (93) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1997, 23, 327. (94) Chaplin, M. http://www.lsbu.ac.uk/water/. (95) Nada, H.; Furukawa, Y. J. Cryst. Growth 2005, 283, 242–256. (96) Case, D. A.; Darden, T. A.; Cheatham, T. E. I.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Merz, K. M.; Wang, B.; Pearlman, D. A.; Crowley, M.; Brozell, S.; Tsui, V.; Gohlke, H.; Mongan, J.; Hornak, V.; Cui, G.; Beroza, P.; Schafmeister, C.; Caldwell, J. W.; Ross, W. R.; Kollman, P. A. AMBER 8.0; University of California, San Francisco: San Francisco, CA, 2004. (97) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33–38. (98) Dang, L. X. J. Chem. Phys. 1992, 97, 2659–2660. (99) Grishina, N.; Buch, V. J. Chem. Phys. 2004, 120, 5217–5225. (100) Kolafa, J.; Onc
ak, M. J. Phys. Chem. C 2010, 114, 20518–20522. (101) Kroes, G. J. Surf. Sci. 1992, 275, 365. (102) Picaud, S. J. Chem. Phys. 2006, 125, 174712. (103) Berendsen, H. J. C.; Grigera, J. R.; Straasma, T. P. J. Phys. Chem. 1987, 91, 6269. (104) Soper, A. K.; Phillips, M. G. Chem. Phys. 1986, 107, 47.

5982

dx.doi.org/10.1021/jp110391q |J. Phys. Chem. A 2011, 115, 5973–5982