1 May 1998
Chemical Physics Letters 287 Ž1998. 429–434
Determination of the structure and stability of water clusters using temperature dependent techniques Ingrid M. Quintana, Wilfredo Ortiz, Gustavo E. Lopez ´
)
Department of Chemistry, UniÕersity of Puerto Rico at Mayaguez, P.O. Box 9019, Mayaguez, Puerto Rico 00681-9019 Received 3 December 1997; in final form 2 February 1998
Abstract The lowest energy equilibrium structure of various water clusters, ŽH 2 O. n , where n s 2–9, was determined using simulated annealing techniques. The interparticle potential was represented by the TIP4P potential and the thermodynamic stability of the systems was computed using a method based on re-scaling the interparticle potential in order to solve the appropriate equations. The structures obtained are similar to the ones reported by other researchers using first principle techniques. All the systems considered are thermodynamically stable in a temperature range where the cluster is in a solid-like form. q 1998 Elsevier Science B.V. All rights reserved.
1. Introduction The study of molecular clusters has been an active area of research over the past years because of their importance in many chemical, physical, biological and environmental processes. Due to their unique behavior, clusters present attractive possibilities for innovative technological applications. Moreover, their detailed characterization helps in the understanding of the behavior of bulk systems. An example of this is the study of water clusters, which have been directly implicated in several problems such as the formation of acid rain, the absorption of sunlight by clouds, and the nucleation of water droplets, to mention some. For these reasons water clusters have been extensively studied from both experimental and theoretical points of view over the past few years. From an ) Corresponding author. Fax: q1-787-265-3849; e-mail:
[email protected]
experimental line of research, the development of high-resolution spectroscopy has allowed the characterization of vibrational and rotational spectra of water clusters. Specifically, Saykally et al.w1–4x have performed far-infrared laser vibration-rotational tunneling experiments on supersonically cooled water clusters. These studies provide information about the structures of various water clusters, such as the trimer, tetramer, pentamer and hexamer. From a theoretical direction, many different models and approaches, such as ab initio techniques, have been applied to the calculation of the unusual properties of condensed water in terms of the properties of its constituent molecules. Results from these studies w5,6x are in good agreement with the ones obtained experimentally. Moreover, detailed characterization of the potential energy surface of small water clusters has been discussed w7x and compared to ab initio results. All of the above studies have provided important insights into the cooperativity effects in hydrogen bonding, aqueous solvation and
0009-2614r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 Ž 9 8 . 0 0 1 6 7 - 5
430
I.M. Quintana et al.r Chemical Physics Letters 287 (1998) 429–434
hydrogen-bond network rearrangement. Previous theoretical studies have also considered the thermodynamic properties of water clusters using temperature dependent formalisms. Tsai and Jordan w8,9x characterized the solid–liquid phase transition for ŽH 2 O. 8 using various Monte Carlo techniques and interparticle potentials. A solid–liquid coexistence was clearly observed for the potentials used and explained in terms of the energies of low-lying stationary points in the potential energy surface. Similarly, Wales and Ohmine w10x studied the coexistence between solid-like and liquid-like clusters using molecular dynamics simulations. Loops in the microcanonical ensemble and bimodalities in the probability distribution at a given temperature were used for the characterization of the coexistence regions. Most of the previous theoretical studies related to cluster stability are based on temperature independent models, hence stability arguments have centered based on energetic contributions. However, our previous work w11x on weakly bound clusters has demonstrated the importance of the entropic contribution in the stability of certain sized clusters even at very low temperature. Therefore, it is the purpose of the present study to examine the stability of water cluster using a thermodynamic model. Specifically, the lowest energy equilibrium structure of different water clusters is determined using simulated annealing techniques. The characterization of the equilibrium structure is followed by the determination of the thermodynamic stability of those structures at finite temperature. The entropic and enthalpic contributions at low temperatures are established through the computation of numerically exact techniques. The remaining of this paper is as follows: In Section 2 we discuss the theoretical models used to obtain the lowest energy equilibrium structures and thermodynamics properties, and in Section 3 we present and discuss the results obtained. Finally, in Section 4 we summarize our findings.
keeps the geometry of the water monomers fixed and represents the charge distribution by four interaction sites w12,13x. This potential has been termed TIP4P and has been used previously in the description of water clusters. Specifically, the model assumes that the water molecule is composed of three charge centers and one Lennard-Jones ŽLJ. center. The total potential energy, Utotal , in this model is given by the following expression w12x: Utotal s Uel q ULJ q Uc .
Ž 1.
Here, Uel is the sum of the Coulomb interactions between all i and j point charge sites, qi q j Uel s Ý , Ž 2. i , j ri j where qi and q j are point charges i and j, respectively, and ri j is the distance between charges i and j. ULJ is the sum of the Lennard-Jones sites and is given by ULJ s
Aa , b
Ý a ,b
ra12, b
y
Ba , b ra6 , b
,
Ž 3.
where Aab and Ba b are Lennard-Jones parameters, and ri j is the distance between a and b LennardJones sites. In this model, the Lennard-Jones sites a and b are located in the oxygen atoms of the water molecule. Finally, a continuos constraining potential w14x, Uc , is used to confine all atoms in a given space volume, n
Uc s
Ý is1
ž
< ri y R cm < RC
20
/
,
Ž 4.
where R cm is the center of mass of the cluster, 1 n R cm s Ž 5. Ýmr. M is1 i i Here, m i is the mass of particle i, ri is the distance between the center of mass and particle i, M is the total mass, and R C is the constraining radius. 2.2. Simulated annealing
2. Theoretical model 2.1. Interparticle potential The interaction between water molecules is described by a pairwise additive potential model which
As stated, the PES of the water clusters has been explored using simulated annealing techniques within the Metropolis Monte Carlo method. The procedure starts with the selection of an initial state where the rigid water molecules are generated in a random
I.M. Quintana et al.r Chemical Physics Letters 287 (1998) 429–434
configuration. The annealing procedure is performed by decreasing the temperature step by step. The temperature increments are given by DT s ŽTmax y Tmin .rk , where Tmax and Tmin are the maximum and the minimum temperatures of the annealing procedure, respectively, and k is the number of annealing temperatures in the simulation. In this work Tmax s 300 K, Tmin s 5 K and k s 15. The procedure is repeated thousands of times in order to guarantee that a minimum in the PES has been reached. 2.3. Thermodynamic properties The change in Gibbs free energy, DGn , for the process Ž H 2 O . ny1 q H 2 O ™ Ž H 2 O . n Ž 6. is calculated using a procedure similar to that used by Zhang et al. w15x and by Rodrıguez et al. w11x. In ´ this techniques the potential is re-scaled in order to compute the necessary configurational integrals. It can be shown that DGn can be computed using the following expression:
½
DGn s yRT ln Ž PrPc . y ln n y b
1
: intra ld l
H0 ²V
5
,
Ž 7. here, R is the gas constant, T is the absolute temperature of the sample, P is the pressure of a standard state, Pc is the constraining factor, n is the number
431
of water molecules in the cluster, l is a scaling factor, and ² Vintra :l is given by the following expression: ² Vintra :l s
Hd r 3 n Vintra eyb DV
. Ž 8. Hd r 3 n ey b DV The calculation of the Gibbs free energy is reduced to the evaluation of a one-dimensional integral in Eq. Ž7. and the computation of ² Vintra :l, which is calculated using classical Monte Carlo techniques for a given value of l. Assuming that the clusters behave ideally in the vapor phase, the change in enthalpy Ž D Hn . and entropy Ž D Sn . for process Ž6. can be calculated by the following relations: D Hn s DUn y RT , Ž 9. D Sn s Ž D Hn y DGn . rT . Ž 10 . Ž . The one-dimensional integral in Eq. 7 was solved using the Gauss–Legendre method with six integration points. On the other hand, Eq. Ž8. was solved using 500 blocks each one composed of 10 5 warm-up moves and 10 5 moves where data was gathered. The step size of the Monte Carlo simulation was adjusted in order to obtain a 50% acceptance ratio. 3. Results and discussion The lowest energy equilibrium structures are presented in Fig. 1 for ŽH 2 O. n from n s 2 to 9. For
Fig. 1. Lowest energy equilibrium structures for the various water clusters, ŽH 2 O. n , where n s 2 to 9. For visual guidance, thin straight lines show hydrogen bonds.
432
I.M. Quintana et al.r Chemical Physics Letters 287 (1998) 429–434
visual guidance, hydrogen bonds have been shown by thin straight lines. It can be observed that the lowest energy equilibrium structure for the systems considered here are similar to the ones obtained from sophisticated ab initio calculations w5,6x. Namely, the clusters from n s 2 to n s 5 form quasi-planar cyclic structures with the formation of hydrogen bonds. Also, for n s 3, 4 and 5 small out-of-plane distortions are observed. In the case of ŽH 2 O. 6 the lowest energy equilibrium structure is a distorted multimember ring. The failure of the TIP4P potential in obtain the lowest energy isomer for this system, which is a prism-like structure, has been previously discussed in detail by Tsai and Jordan w16x. The clusters from n s 7 to n s 9 show structures that are form of various water rings bound together by hydrogen bonds. In particular, it can be noticed that two four-member ring clusters bound together by hydrogen bonds form the ŽH 2 O. 8 cluster. This cluster has an S 4 symmetry. In order to determine if the systems are thermodynamically stable, we have computed changes in the Gibbs free energy, as well as changes in enthalpy and entropy for reaction Ž6.. Fig. 2 shows the variation in DGn as a function of the number of water molecules in the cluster. As stated, this quantity is associated to the stability of the n-molecule cluster. The various curves represent the thermodynamic quantity as a function of temperature. For the temperatures considered here the clusters remain in a
Fig. 2. Variation of DG as a function of the number of molecular species in the cluster at various temperatures. The data points are connected by straight lines for clarity only. Here, the error bars are within the size of the data points.
Fig. 3. Same as Fig. 2, but the variation of DH as a function of the number of molecular species in the cluster at various temperatures.
solid-like structure. In all cases and in the whole temperature range studied, DGn is negative which implies that all the clusters are thermodynamically stable. Moreover, two minimum in DGn for the n s 4 and n s 8 clusters can be observed. Hence, magic number for water clusters can be identified when four-member quasi-planar rings are formed and bound by hydrogen bonds. It is also observed from this figure that as the temperature increases the magic numbers minimum becomes less pronounce, in particular for n s 8. A possible explanation for this is that as the temperature is increased the hydrogen bonds between the two rings in the octamer are broken and isomers associated to two four-member rings not bound by four hydrogen bonds are formed. The stability of these possible isomers should be very similar to the tetramer, a fact that explains the similarity in the Gibbs free energy for both the tetramer and octamer, particularly at a temperature of 160 K. In order to discover the nature of the stabilization of these complexes, the entropic and enthalpic contributions are computed for process Ž6. using Eq. Ž9. and Eq. Ž10., respectively. Fig. 3 shows the variation in D H as a function of the number of water molecules in the clusters. Clearly, a tendency similar to what is observed for the variation in the Gibbs free energy can be identified. Namely, two pronounced minimum can be identified for the tetramer and for the octamer, with the one for the octamer more pronounce. This result clearly implies that the stability
I.M. Quintana et al.r Chemical Physics Letters 287 (1998) 429–434
Fig. 4. Same as Fig. 2, but the variation of DS as a function of the number of molecular species in the cluster at various temperatures.
of these systems is energetically driven as it might be expected owing the fact that hydrogen bonds are being formed. Small variations in D H are observed with variation in temperature. Fig. 4 shows the variation in D S as a function of cluster size at various temperatures. Two main tendencies can be observed from these curves. The first one is for clusters of n - 6 which exhibit negative changes in the entropy and a fluctuating behavior of D S as n increases. Moreover, a small dependence in the values of D S with temperature is observed. On the other hand, for clusters with n ) 6 an increment in the value of D S is observed as the number of molecules forming the clusters increases. A possible explanation for this is that for clusters with n ) 6 the number of possible configurations or isomers increases considerably as the number of molecules increases, hence causing an increment in the change in entropy. However, for the systems with n - 6, which present quasi-planar cyclic structures, the number of configurations is much smaller and hence the entropy does not change much. In most of the cases the values for D S are negative which implies that the clusters are not entropically stabilized. An important feature of the data that can be observed in this figure is the behavior of D S for the various temperatures. Again, for the region where n - 6, small variations in D S are observed as the temperature of the systems is increased. However, for n ) 6 the behavior is quite different, because as the temperature of the cluster is decreased, the value
433
for D S increases. In particular, in the case of n s 9 at T s 80 K the value for D S is positive. A possible explanation for this phenomenon is the following. At high temperatures the thermal motion of the clusters is considerable, causing the similarity in the structures of the various clusters. For example, at high temperatures the n s 8 and n s 9 clusters can be very similar in structure because of the weakly bound nature of the systems and therefore the entropy will be approximately the same Ž D S s 0.. However, at low temperatures, the clusters lost its thermal motion and are more solid-like in nature. When the n s 8 cluster is compared with the n s 9 cluster, both in a more solid-like form, a loss of symmetry can be observed, hence causing an increment in entropy. This result is similar to what was observed for N2 clusters w11x.
4. Conclusions In this study we have found that water cluster systems form ring-like structures for n F 5 and multi-ring structures for larger clusters. The structures obtained are in agreement with previous studies using a variety of theoretical techniques. All the systems considered are thermodynamically stable in the temperature range of T s 80 K to 160 K. The stability is mainly due to an energetic contribution, which is associated to the formation of hydrogen bonds. In particular, magic number clusters where identified when four-member rings are formed. At present we are considering the influence of a flexiblerpolarizable potential in the thermodynamic properties of water clusters. Quantum effects based on a temperature dependent formalism are also being considered.
Acknowledgements A Cottrell College Science Award of the Research Corporation and the NSF-EPSCoR program supported this work. GEL is a Camille and Henry Dreyfus Teacher Scholar.
434
I.M. Quintana et al.r Chemical Physics Letters 287 (1998) 429–434
References w1x K. Liu, J.D. Cruzan, R.J. Saykally, Science 271 Ž1996. 929. w2x N. Pugliano, R.J. Saykally, Science 257 Ž1992. 1937. w3x J.D. Cruzan, L.B. Braly, K. Liu, M.G. Brown, J.G. Loeser, R.J. Saykally, Science 271 Ž1996. 59. w4x K. Liu, M.G. Brown, J.D. Cruzan, R.J. Saykally, Science 271 Ž1996. 929. w5x S.S. Xantheas, T.H. Dunning, J. Chem. Phys. 99 Ž1993. 8774. w6x C. Lee, H. Chen, G. Fitzgerald, J. Chem. Phys. 102 Ž1995. 1266. w7x C.J. Tsai, K.D. Jordan, J. Chem. Phys. 97 Ž1993. 11227.
w8x w9x w10x w11x w12x w13x w14x w15x w16x
C.J. Tsai, K.D. Jordan, J. Chem. Phys. 95 Ž1991. 3850. C.J. Tsai, K.D. Jordan, J. Chem. Phys. 99 Ž1993. 6957. D.J. Wales, I. Ohmine, J. Chem. Phys. 98 Ž1993. 7245. I. Rodrıguez, A.J. Acevedo, G.E. Lopez, Mol. Phys. 90 ´ ´ Ž1997. 943. L. Perera, M.L. Berkowitz, J. Chem. Phys. 95 Ž1991. 1954. W.L. Jorgensen, J.D. Madura, Mol. Phys. 56 Ž1985. 1381. D.L. Freeman, J.D. Doll, Adv. Chem. Phys. B 70 Ž1988. 139. C. Zhang, D.L. Freeman, J.D. Doll, J. Chem. Phys. 91 Ž1989. 2489. C.J. Tsai, K.D. Jordan, Chem. Phys. Lett. 213 Ž1993. 181.