THE JOURNAL OF CHEMICAL PHYSICS 125, 094508 共2006兲

Three-dimensional square water in the presence of an external electric field M. Girardia兲 Instituto de Física, Universidade de São Paulo, C.P. 66318, 05315-970 São Paulo, São Paulo, Brazil

W. Figueiredob兲 Departamento de Física, Universidade Federal de Santa Catarina, 88040-900 Florianópolis, Santa Catarina, Brazil

共Received 13 June 2006; accepted 11 August 2006; published online 6 September 2006兲 In this work we study a tridimensional statistical model for the hydrogen-bond 共HB兲 network formed in liquid water in the presence of an external electric field. This model is analogous to the so-called square water, whose ground state gives a good estimate for the residual entropy of the ice. In our case, each water molecule occupies one site of a cubic lattice, and no hole is allowed. The hydrogen atoms of water molecules are disposed at the lines connecting nearest-neighbor sites, in a way that each water can be found in 15 different states. We say that there is a hydrogen bond between two neighboring molecules when only one hydrogen is in the line connecting both molecules. Through Monte Carlo simulations with Metropolis and entropic sampling algorithms, and by exact calculations for small lattices, we determined the dependence of the number of molecules aligned to the field and the number of hydrogen bonds per molecule as a function of temperature and the intensity of the external field. The results for both approaches showed that, different of the two-dimensional case, there is no maximum in the number of HBs as a function of the electric field. However, we observed nonmonotonic behaviors as a function of the temperature of the quantities of interest. We also found the dependence of the entropy on the external electric field at very low temperatures. In this case, the entropy vanishes for the value of the external field for which the contributions to the total energy coming from the HBs and the field become the same. © 2006 American Institute of Physics. 关DOI: 10.1063/1.2348866兴 I. INTRODUCTION

It is well known that the hydrogen-bond 共HB兲 network present in liquid water plays an important role in its thermodynamical and also dynamical properties.1,2 High dielectric constant, negative change of volume on melting, and high melting and boiling points are some examples of the unusual behavior consequence of hydrogen bonding. The reactivity, viscosity, and diffusibility are also affected, lowering with the strengthening of the HBs. By imposing an external electric field one can destroy the HB network and change considerably the properties of water. Some recent papers focused attention on this fact and studied the influence of the electric field on the dielectric constant, polarization, and the crystallization of supercooled water to proton ordered ice forms, known as electrofreezing.3–10 Sutmann employed molecular dynamics simulations on extended simple point charge 共SPC/E兲 water under strong electric fields and obtained a completely polarized icelike structure at room temperatures. This result was confirmed by Yeh and Berkowitz, which also calculated the dielectric constant as a function of the external electric field E, showing that it decreases as 1 / E. Other dynamical and structural properties were obtained by Choi et al.11 and by Vegiri and Schevkunov6 and Vegiri7 for bulk and small water clusters under strong electric fields. Using a density funca兲

Electronic mail: [email protected] Electronic mail: [email protected]

b兲

0021-9606/2006/125共9兲/094508/6/$23.00

tional theory, Choi et al. observed that, depending on the relative direction of the applied field, the HBs are weakened or strengthened, for the case of cyclic and linear water clusters containing three to five molecules. Electric fields perpendicular to the plane of the rings of water clusters increase the HB length, decreasing their strength, although fields parallel to the linear clusters enhance the hydrogen bonding due to the reduction of their average length. By employing molecular dynamics for a TIP4P water model, Vegiri7 determined the structural and reorientational relaxation times for liquid water and additionally found the self-diffusion coefficients as a function of the external electric field. An increasing relaxation time and decreasing diffusion coefficient with the field were found, together with an spatial anisotropy, where the diffusion in the plane perpendicular to the field is enhanced if compared to that along the field direction. In the present work we extended to three dimensions an earlier studied two-dimensional model12–15 for the hydrogenbond network under an external electric field. Here, a cubic lattice is completely filled with water molecules, each one occupying a single site, and assuming different orientations 共states兲, forming up to four hydrogen bonds with their nearest neighbors. Applying Monte Carlo simulations and exact calculations for small lattices, we obtained the number of HBs per molecule and the fraction of aligned molecules to the field as a function of temperature and magnitude of the external field. The residual entropy and the size distribution of clusters of molecules aligned to the field were also calcu-

125, 094508-1

© 2006 American Institute of Physics

Downloaded 19 Aug 2008 to 200.18.45.248. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

094508-2

J. Chem. Phys. 125, 094508 共2006兲

M. Girardi and W. Figueiredo

FIG. 1. The 15 different states of the bonding solvent particles on a cubic lattice. Note that two neighboring molecules form a bond if there is only one arrow between them, as for the pair of particles in states 1 and 2 shown in the figure.

lated. In this model, the direction of the external field was chosen to favor a given orientation of the water electric dipoles, and depending on this choice, the overall behavior of the system seems to be completely different. A comparison between the results in two and three dimensions is also made, including a discussion about a possible phase transition induced by the external field at low temperatures. II. MODEL AND CALCULATIONS

The square water model extended to three dimensions is represented by a cubic lattice of linear size L where each site is occupied by a water molecule, one at each site. No hole is allowed, and consequently no density anomaly is present. A water molecule can be in one of its 15 different states representing the possible dispositions of the two hydrogen atoms in the six directions of the lattice. All possible states are shown in Fig. 1. Water molecules interact via directional HB-like coupling, where two neighboring water molecules form a HB if there is only one hydrogen in the line connecting both molecules. We associate the energy eHB = −␧ 共␧ ⬎ 0兲 for each HB. The molecules also interact with an external electric field that favors energetically one given state of water. In this way, each aligned molecule contributes with eh = −h to the total energy, which can thus be written as U = −␧NHB − hNh, where NHB is the total number of hydrogen bonds and Nh is the total number of molecules aligned to the external field. Similarly to the two-dimensional version of the model, for a vanishing value of the field and nonzero temperature, there is no phase transition, distortion, or density fluctuations, and, for all temperatures, the hydrogenbond network percolates. The number of HBs varies from 2 per molecule at zero temperature to 4 / 3 per molecule at very high temperatures. Since for T → ⬁ the molecules are uncorrelated, the probabil5 共1 − 155 兲, ity of two neighbors to form a HB is PHB = 2 ⫻ 15 5 where the factor 15 represents the fraction of states which has one arrow pointing to a given neighbor, and the factor 2 indicates that the arrow can point in both directions. Then, the mean number of HBs per molecule is nHB = 3PHB = 4 / 3, since the bond can be on each one of the directions x, y, and z. Two different approaches were employed to study the thermodynamics of the system. The first one is the exact

calculation of the partition function for a finite small system. In this case, we generate all configurations for a lattice with linear size L = 2 and periodic boundary conditions. We write down the partition function as Z = 兺 Nab exp关␤共a␧ + bh兲兴, a,b

where Nab is the number of states with a hydrogen bonds and b molecules aligned to the external field, ␤ = 共kBT兲−1, kB is the Boltzmann constant, and T is the temperature. The temperature T is measured in units of ␧ / kB. For L = 2, the indices in the above sum are in the range of a = 0 , . . . , 16 and b = 0 , . . . , 8. Since there is no critical phenomenon in this system, finite-size effects are not relevant, and even small lattices can predict well its thermodynamics. This fact will become evident when we compare the exact results for L = 2 with the simulation data for bigger lattice sizes. The ¯ , mean numquantities of interest as the mean total energy U ber of HBs NHB, mean number of aligned molecules Nh, and the entropy are given by ¯ = − ⳵ ln共Z兲 , U ⳵␤ Nh = −

NHB = −

1 ⳵ ln共Z兲 , ␤ ⳵h

1 ⳵ ln共Z兲 , ␤ ⳵␧

S = − ␤2

⳵ ln共Z兲 . ⳵␤ ␤

Another method we used to study the present model is the Monte Carlo method. Here we performed simulations with the Metropolis16 and entropic sampling17 algorithms. In both cases, a lattice with linear size L = 16 is initially filled with water molecules in random states. Then, we allow the system to evolve in time by changing the states of each molecule. In one Monte Carlo step 共MCS兲, we try to change the states of all molecules. If the system evolves via Metropolis algorithm, the acceptance rate for a change from state ␯i to ␯ j is A共␯i → ␯ j兲 = min关1 , exp共−␤⌬U兲兴, where ⌬U is the change in energy. After typically 104 MCS the system reaches the thermal equilibrium and the quantities of interest can be calculated. In order to obtain the entropy of the system, we also used the entropic sampling algorithm. In this method, the changes are accepted with probability,

Downloaded 19 Aug 2008 to 200.18.45.248. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

094508-3

Water in the presence of an electric field

冦

if ␳n共U␯i兲 ⱖ ␳n共U␯ j兲

1

A共␯i → ␯ j兲 = ␳n共U␯i兲

␳n共U␯ j兲

otherwise,

冧

where ␳n is the density of states 共DOS兲 in the nth selfconsistent step of the simulation. This DOS is updated as follows: at the beginning of the simulation, ␳ is defined as being flat, so ␳0共U兲 = 1 for all U, and the acceptance rate is equal to 1, until the next update of ␳. While the system evolves in time, initially performing a random walk in the phase space, we build a temporary histogram of energies g共U兲 updated at each MCS. After g共U兲 has accumulated, in average, at least 100 units for each bin in which the energy was divided 共this number is arbitrary and must be chosen by the experience兲, we start to update ␳ as ln关␳n+1共U兲兴 =

再

ln关␳n共U兲g共U兲兴 if g共U兲 ⬎ 0 ln关␳n共U兲兴

if g共U兲 = 0,

J. Chem. Phys. 125, 094508 共2006兲

冎

for each value of U. For n ⬃ 106, ␳ converges, and the system performs a random walk in the space of energies. As a characteristic of the entropic sampling algorithm, the distribution ␳共U兲 is related to the entropy S共U兲 by the relation ␳共U兲 ⬀ exp关S共U兲兴. This means that the knowledge of the entropy for some value of the energy is sufficient to determine the former for any other arbitrary value of energy 关and also for any temperature, since ␤ = ⳵US共U兲兴. In our case, the entropy for T → ⬁ is simply L3 ln 15, since each molecule can be in one of the 15 states. The calculation of S共T兲 is then straightforward. III. RESULTS AND DISCUSSIONS

First, we will consider an external electric field that favors a given state in the range of 1 and 12 共here we will choose state 1兲. In this case, a pair of neighboring molecules in a given state of this range can form an HB. The same does not occur for a pair of molecules in one of the states 13, 14, and 15 that never forms an HB no matter their relative positions. Figure 2 exhibits the number of hydrogen bonds per molecule and the fraction of aligned molecules to field as a function of the temperature and the external field when state 1 is the favored one. As in the two-dimensional square water,12 nHB = NHB / N decreases with the temperature 关Fig. 2共a兲兴 for any value of h, and the same happens with the fraction of aligned molecules nh. At nonzero temperatures, the external field increases the number of molecules in state 1 关Fig. 2共d兲兴 and consequently the number of HBs. This occurs because a system with all molecules in state 1 is fully bounded 共nHB = 2兲. The same behavior was obtained by Kiselev and Heinzinger for the SPC/E water under applied electric field,18 and by Sutmann10 for the Bopp-JancsóHeinzinger model19 of water in the presence of a strong electric field. In the work of Sutmann,10 the increase in the number of four-bounded molecules leads to a freezing of the translational and rotational degrees of freedom and a structuring of the lattice. Suresh et al.8 developed a theoretical framework for water in strong electric fields, which also predicts an increase in the number of HBs with the field. The

FIG. 2. Number of hydrogen bonds per molecule nHB and fraction of aligned molecules as a function of temperature 关共a兲 and 共c兲兴 and external field 关共b兲 and 共d兲兴. Here state 1 is the favored one. 共a兲 and 共c兲 are simulations for L = 16 and h = 0 共circles兲, h = 4 共triangles兲, and h = 6 共crosses兲. Dashed lines are the exact results for L = 2. 共b兲 and 共d兲 are simulations for L = 16 共circles兲 and exact results for L = 2 共dashed line兲 for T = 1.

authors justify this behavior by the reduction of the thermal motion and the enhancement of the alignment favoring the bonding. Also, note in Fig. 2 the good agreement between simulation data for L = 16 and the exact results for L = 2. As a matter of comparison, we exhibit in Fig. 3 the results for the number of HBs per molecule at zero external field obtained by Suresh and Naik20 from the experimental values of the dielectric constant of water and our simulation data. As stated by Nadler and Krausche,13,14 a direct comparison of the square water model with real water is only possible by rescaling the temperature. In this way, the hydrogen-bond energy ␧ must be viewed as an effective free energy ␧ = ␧⬘ − T␴, where ␧⬘ is the estimated hydrogen-bond energy and ␴ is an entropy factor. This generalization is necessary in order to account for entropy differences between open and closed HBs due to not considered degrees of freedom 共e.g., vibrational兲. Here we used ␧⬘ = 1.5⫻ 10−20 J and

FIG. 3. Number of hydrogen bonds per molecule nHB as a function of temperature in kelvin. Simulations for L = 16, h = 0 and rescaled ␧ to ␧ = 1.5⫻ 10−20 – 3.3⫻ 10−23T 共dashed line兲 and the results of Suresh and Naik 共continuous line兲.

Downloaded 19 Aug 2008 to 200.18.45.248. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

094508-4

M. Girardi and W. Figueiredo

FIG. 4. Number of hydrogen bonds per molecule nHB as a function of the external electric field E given in V/m. Exact results for L = 2, T = 289 K 共connected circles兲, T = 350 K 共connected triangles兲 and rescaled ␧ to ␧ = 1.5⫻ 10−20 – 3.3⫻ 10−23T. The results of Suresh et al. for T = 289 K 共continuous line兲 and T = 350 K 共dashed line兲.

␴ = 3.3⫻ 10−23 J / K.13,14 For these values of the parameters, our results are in accordance with those reported on the work of Suresh and Naik.20 In the case of fixed temperature and variable external electric field, we compare our results with those obtained by Suresh et al.8 for a theoretical model of water. The mean HB number for two different temperatures are shown in Fig. 4, where we used a constant value for the magnitude of the water dipole moment, ␮ = 1.854 D. Here we also defined E = h / ␮. As can be seen in Fig. 4, for both studied temperatures, the behavior is only qualitatively similar to that obtained by Suresh et al. In our case, water is completely polarized at lower fields and is more susceptible to them. This can be explained by the limited translational and rotational degrees of freedom present in our model, which made it easier to align the molecules. A different scenery arises when the external field favors the states such as 13, 14, and 15. Since a system with all molecules in one of these states does not present hydrogen bonds, a nonmonotonic behavior of nh and nHB appears. In Fig. 5 we show that the fraction of HBs decreases with the temperature for low values of h and increases for high values. At intermediate values of the external field there are local maxima and minima in the curves nHB vs T. This change of behavior occurs at h = 4, for which the energy of a nonbounded aligned molecule to the field is the same as the one forming four HBs. Also in Fig. 5, we exhibit the exact results for L = 2 for the same values of the external field. Note that for h = 0 and h = 6, the exact results for the small lattice is quantitatively similar to those obtained by Monte Carlo simulations. However, at values near h = 4, the agreement between exact and simulation results is only qualitative. In this case, the ordered state does not fit very well to a small lattice, and finite-size effects become important. At h = 4 and at zero temperature, the system suffers a first-order phase transition, in which we observe a coexistence between a phase where all molecules are in state 15 and

J. Chem. Phys. 125, 094508 共2006兲

FIG. 5. Number of hydrogen bonds per molecule nHB as a function of temperature. Simulations for L = 16 and h = 0 共circles兲, h = 3.9 共squares兲, h = 4.0 共triangles兲, h = 4.1 共crosses兲, and h = 6 共stars兲. The dashed lines are the exact results for L = 2. Here the external field favors state 15.

another one, where only half of the molecules are in this state. At T = 0, the free energy of the system is exactly its total energy, which is given by U / N = −2␧ − h / 2 for 0 ⬍ h ⬍ 4, U / N ⯝ −1.6␧ − 0.6h for h = 4, and U / N = −h for h ⬎ 4. The first derivative of U is not continuous at h = 4, and a first-order phase transition takes place. The coexistence between the two phases can be appreciated in Fig. 6 that gives the size distribution of clusters of aligned molecules. For this choice of parameters, a large aggregate containing almost all molecules in state 15, which we chose as the favored one, is surrounded by isolated ones. The fraction of aligned molecules is nh ⯝ 0.6. Another first-order phase transition also occurs at T = 0 and h = 0. At h = 0−, nh = 0, while at h = 0+, nh = 1 / 2. Exactly at h = 0 there is a coexistence of phases with the fraction of aligned molecules equal to 1 / 15. Again, the total energy has a cusp at this point, changing its slope from zero for h 艋 0 to −1 / 2 for 0 ⬍ h ⬍ 4, which signals a firstorder phase transition. The behavior at T = 0 can be summarized as follows:

FIG. 6. Size distribution curve for clusters of molecules aligned to the field. Simulation results for T = 0, L = 16, where the state 15 is the favored one.

Downloaded 19 Aug 2008 to 200.18.45.248. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

094508-5

J. Chem. Phys. 125, 094508 共2006兲

Water in the presence of an electric field

FIG. 7. Fraction of aligned molecules nh and number of HBs per molecule nHB as a function of the external field for T = 0.2. Simulation 共circles for nh and crosses for nHB兲 and exact results for L = 2 共dashed lines兲. Here state 15 is the favored one.

nHB = 2 and nh = 1 / 15 for h = 0, nHB = 2 and nh = 1 / 2 for 0 ⬍ h ⬍ 4, nHB ⯝ 1.6 and nh ⯝ 0.6 for h = 4, and nHB = 0 and nh = 1 for h ⬎ 4. For nonzero temperatures, as we can see in Fig. 7, nHB is a monotonic decreasing function of the field magnitude, and there is no local maximum for some value of h, different from the two-dimensional version of the model where this maximum is present for a wide range of temperatures. The weakening of the HBs with the external field was observed in molecular dynamics simulations for small linear water clusters at fields perpendicular to the cluster.11 There, the electric field polarizes the water molecules in the direction perpendicular to the HBs, decreasing the dipole-dipole interaction and, consequently, the hydrogen bonds. In Fig. 8 we plot the fraction of aligned molecules as a function of temperature, and again a nonmonotonic behavior arises. For example, at h = 3.9, the competition between the

FIG. 9. Entropy per particle S / L3 as a function of the external field for T = 0.2. Simulations for L = 6 共circles兲 and exact results for L = 2 共continuous line兲. Here state 15 is the favored one. Inset: linear behavior of the total entropy vs L3 共circles兲 for h = 2.

energies to align a molecule to the field and that to form a hydrogen bond is clear. At T = 0, half molecules are aligned and nHB = 2. As we increase the temperature, nHB must decrease and, at sufficiently high values of h, the best way to do it is to align bounded molecules, breaking their bonds. Finally, in Fig. 9, we present the simulation data for the entropy per molecule at T = 0 as a function of the external electric field. Note that there is a residual entropy for h 艋 4, which implies in a highly degenerated energy minimum. In two dimensions, the residual entropy is nonzero only at h = 0 and h = 4. The inset of Fig. 9 displays the total entropy as a function of the volume of the system for h = 2, whose linear behavior indicates a nonvanishing entropy in the thermodynamic limit. In conclusion, we have performed Monte Carlo simulations and exact calculations for a three-dimensional model of the water hydrogen-bond network under an external electric field. We obtained the fraction of hydrogen bonds and aligned molecules as a function of temperature and the magnitude of the external field, for two different favored orientations of the water molecules. Depending on the favored state, a decreasing or increasing number of HBs with the field is observed, in accordance with molecular dynamics simulations and density functional calculations for small water clusters under electric fields applied in different directions. We also calculated the residual entropy, which is nonzero for h 艋 4, and found two first-order phase transitions at T = 0. The coexisting phases are characterized by different densities of aligned molecules. ACKNOWLEDGMENTS

We would like to thank the financial support of the Brazilian agencies FAPESP, CNPq, and FAPESC. F. Franks, Water: a Comprehensive Treatise 共Plenum, New York, 1973兲. D. Chandler, Nature 共London兲 417, 491 共2002兲. 3 W. Sun, Z. Chen, and S. Huang, Fluid Phase Equilib. 238, 50 共2005兲. 4 S. Harinipriya and M. V. Sangaranarayanan, Indian J. Chem. Section A 42, 965 共2003兲. 1

FIG. 8. Fraction of aligned molecules nh as a function of temperature. Simulations for L = 16 and h = 0 共circles兲, h = 3.9 共squares兲, h = 4.0 共triangles兲, h = 4.1 共crosses兲, and h = 6 共stars兲. The dashed lines are the exact results for L = 2. Here the state 15 is the favored one.

2

Downloaded 19 Aug 2008 to 200.18.45.248. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp

094508-6

J. Chem. Phys. 125, 094508 共2006兲

M. Girardi and W. Figueiredo

I. Yeh and M. L. Berkowitz, J. Chem. Phys. 110, 7935 共1999兲. A. Vegiri and S. Schevkunov, J. Chem. Phys. 115, 4175 共2001兲. 7 A. Vegiri, J. Mol. Liq. 112, 107 共2004兲. 8 S. J. Suresh, A. V. Satish, and A. Choudhary, J. Chem. Phys. 124, 074506 共2006兲. 9 S. J. Suresh, J. Chem. Phys. 122, 134502 共2005兲. 10 G. Sutmann, J. Electroanal. Chem. 450, 289 共1998兲. 11 Y. C. Choi, C. Pak, and K. S. Kim, J. Chem. Phys. 124, 094308 共2006兲. 12 M. Girardi and W. Figueiredo, J. Chem. Phys. 117, 8926 共2002兲.

W. Nadler and T. Krausche, Phys. Rev. A 44, R7888 共1991兲. T. Krausche and W. Nadler, Z. Phys. B: Condens. Matter 86, 433 共1992兲. 15 E. H. Lieb, Phys. Rev. 162, 162 共1967兲. 16 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 共1953兲. 17 J. Y. Lee, Phys. Rev. Lett. 71, 2353 共1993兲. 18 M. Kiselev and K. Heinzinger, J. Chem. Phys. 105, 650 共1996兲. 19 P. A. Bopp, G. Jancsó, and K. Heinzinger, Commun. Pure Appl. Math. 98, 129 共1983兲. 20 S. J. Suresh and V. M. Naik, J. Chem. Phys. 113, 9727 共2000兲.

5

13

6

14

Downloaded 19 Aug 2008 to 200.18.45.248. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp