Monte Carlo Investigation of the Thermodynamic Properties of (H2O)n

ARTICLE pubs.acs.org/JPCB

Monte Carlo Investigation of the Thermodynamic Properties of (H2O)n and (H2O)nH2 (n = 220) Clusters Glen L. Holden*,† and David L. Freeman Department of Chemistry, University of Rhode Island, Kingston, Rhode Island 02881, United States ABSTRACT: The classical thermodynamic properties of water clusters modeled by the TIP4P potential with and without the inclusion of molecular hydrogen are calculated using Monte Carlo methods. The temperature-dependent heat capacity curves of the dimer and the smaller pure water clusters having ring structures show low-temperature anomalies arising from the onset of a transition from librational motion to free rotational motion. Pure water clusters having cage structures display heat capacity anomalies characteristic of “melting” phase changes. The addition of molecular hydrogen to the water clusters has little impact on the structure of the water core, but low-temperature heat capacity anomalies are observed that arise from the onset of hindered to free rotational and translational motion of the hydrogen molecule with respect to the skeletal moieties of the core. The Gibbs free energy changes associated with the growth of pure water clusters and the addition of molecular hydrogen to the water clusters are determined using a combination of state-integration methods and an application of the GibbsHelmholtz equation. For the temperature integration of the GibbsHelmholtz equation, a quadrature method is introduced that avoids numerical difficulties arising from singularities in the integrand at low temperatures. For the growth of pure water clusters, the fine structure of the enthaply and the low-temperature Gibbs free energy as a function of cluster size is rationalized using an ansatz of molecular dipole orientations. For the addition of molecular hydrogen to water clusters, the Gibbs free energy change is a virtually flat function of cluster size showing no fine structure.

1. INTRODUCTION There has been continuing interest116 in the properties and behaviors of clusters of water molecules. As with other finite systems, water clusters are both inherently interesting and have been examined as models for the behavior of bulk water. Because water comprises approximately 70% of the earth’s surface and is essential for life, water is arguably the most important of all chemical substances. From a scientific perspective, water has intriguing properties that have inspired many experimental and theoretical investigations. Examples of the special properties of water include the low density of the solid phase compared to the liquid phase, the complex phase diagram of the solid state, the high dielectric constant of water, and the residual third-law entropy. Many of the special properties of water can be understood from the unique strength of the hydrogen bonding found in the condensed phase along with the substantial dipole moments of the individual water molecules. The special properties of water are well-known and partially understood in the bulk phase. It is of interest to examine the extent to which analogues of some of these bulk anomalies are also observed in water clusters. Water clusters are also important owing to their role in condensation phenomena. As with any first-order phase transition, the dynamics of water condensation is governed by nucleation. In classical nucleation theory,17 the free energy of formation of the clusters that are the intermediate nuclei in the nucleation process provide important information concerning nucleation r 2011 American Chemical Society

rates. Understanding the thermodynamic properties of the cluster intermediates is essential to develop a full picture of water condensation. In the current work, we study computationally the thermodynamic properties of water clusters ranging in size from n = 2 water molecules to 20. Of particular interest to us are the phase change properties and the Gibbs free energies of the cluster growth process. Our vehicle for studying the phase change properties is the heat capacity of the clusters. For some cluster sizes, the cluster heat capacity as a function of temperature has been investigated previously5,7,8,10,13,14,16 using a variety of potential models. These studies have shown a strong sensitivity of the location of heat capacity peaks as a function of temperature to the nature of the potential model used.5 Perhaps the most accurate approach to the choice of potential model would be to use ab initio calculations and determine the energy of the system on the fly. Unfortunately, such an approach is prohibitively expensive for heat capacity determinations especially for the larger clusters. Owing to quasiergodicity problems, statistically accurate heat capacity calculations can require tens of millions of Monte Carlo energy evaluations, something that is not currently possible using available computer technology with on-the-fly ab Received: February 1, 2011 Revised: March 4, 2011 Published: April 04, 2011 4725

dx.doi.org/10.1021/jp201082p | J. Phys. Chem. B 2011, 115, 4725–4744

The Journal of Physical Chemistry B initio evaluation of the energy. In this work we use the TIP4P potential18 which is believed to be satisfactory for clusters smaller than n = 17.19 Other, more sophisticated, potential models have been used in previous work, but only for limited sets of small clusters.20 We associate heat capacity peaks for the larger clusters with solid to liquidlike phase changes. Unlike some other studies,8,21 we make no attempt to identify the phase changes using the Landau free energy, and we examine the heat capacity anomalies as reflections of the structure of the underlying potential energy surface. Because of our interest in the relation between the potential surface and the observed features of the heat capacity curves and other thermodynamic properties, we extend our investigations to lower temperatures than previous studies. While the low temperatures studied here are in regions that physically may require a quantum treatment, we find intriguing information about the potential surface by examining the classical heat capacities at these low temperatures. Quantum effects in some other systems can completely remove some phase change behaviors. For example, in Ne38 not only do quantum effects remove the interesting solidsolid phase change, but the entire melting peak in the heat capacity observed in classical simulations is absent when quantum effects are included.22 Similar dramatic changes in thermodynamic properties with the inclusion of quantum effects have been observed22 in the SPC model of (H2O)8. Nonetheless, the classical studies remain interesting, because the classical features reveal much about the underlying potential energy surface. To calculate the Gibbs free energies, we use a modified stateintegration method that combines the Kirkwood method23 with integration of the GibbsHelmholtz equation. The combination is useful, because the Kirkwood scaling method requires a large quadrature grid resulting in significant computational overhead. By generating a specialized and novel quadrature rule, we are able to perform the GibbsHelmholtz integration with modest computational cost. At low temperatures, we find the Gibbs free energy as a function of cluster size has interesting fine structure. By resolving the Gibbs free energy into its enthalpic and entropic components, we are able to generate a useful model for the energetics of the clusters that includes the observed orientations of the dipoles. Mixed water clusters have been found to be useful models for aqueous solvation.2427 Because solvated ions play an essential role in chemistry, and because charges simplify experimental observations, much of the work on mixed water clusters has focused on solvated cations and anions. Much of the structural work on neutral clusters of single-molecule solutes in water clusters has been motivated by using the clusters to understand better the properties of clathrate hydrates. Examples include studies of methane solvated in water cages including the effect of pressure.28 Motivated by experimental and theoretical studies of hydrogen clathrate hydrates,29 in addition to studying the thermodynamic properties of pure water clusters, we examine (H2O)nH2. As discussed later for both pure water and molecular hydrogen with water clusters, we include an external constraining potential about the center of mass. Such constraining potentials, included to handle evaporation events in simulations, do induce an internal pressure on the cluster.31 Unlike some of the available work on methane in water clusters,28 we do not vary this internal pressure, and we choose constraining radii having minimal effects on the thermodynamic properties while enabling ergodic simulations.

ARTICLE

Table 1. Partial Charges and Lennard-Jones Interaction Parameters for TIP4P Water and the Rigid H2 Molecule Used in MC Simulationsa atom

q (e)

σii (bohr)

εii (kJ mol1) 0.6490

O(H2O)

0.0000

5.9594

H(H2O)

þ0.5200

0.0000

0.0000

M(H2O)

1.040

0.0000

0.0000

H(H2)

þ0.4932

0.0000

0.0000

Hcm(H2)

0.9864

5.74099

0.2852

a

The site labeled M represents the displaced charge 0.15 Å from the oxygen atom of the water molecule.

As we find later, with such constraining potentials, the effect of the hydrogen molecule on the structure of the water cage is minimal. However, we discover intriguing physical effects from the motion of the hydrogen molecule on the low-temperature thermodynamic features. The contents of the remainder of this paper are as follows. In the next section we discuss the computational methods used and the potential models chosen for the waterwater and waterhydrogen interactions. The results are presented in section 3, and in the final section we discuss our findings.

2. METHODS Monte Carlo methods represent the principle computational tool used in the current work. Because the Monte Carlo approaches used are based on the Metropolis method (MMC)32 and are unexceptional, we make no effort to review the principal techniques here. Owing to well-known ergodicity problems in simulating the properties of clusters,33,34 especially in the “melting region” of temperatures, we also use parallel tempering methods3537 that have been shown to be highly useful in other cluster studies.33 Details concerning the parameters used in the current work are discussed in section 3. In the remainder of this section we outline those aspects of the methods used that provide the computational details needed to reproduce the work. 2.1. Model Potential. The water clusters in the present work are simulated using the classical TIP4P water potential.18 A structural comparison between water clusters, over a range of sizes, between the simpler TIP4P and the physically richer TTM2-F potential has been reported by Hartke.19 Comparisons between the structures of TIP4P water clusters and other models show good accuracy for TIP4P clusters, at least for sizes n e 16.19 In the TIP4P model, a negative charge is placed 0.15 Å from the oxygen atom, along the C2 axis bisecting the HOH angle of the water molecule. The parameters for the TIP4P water potential are listed in Table 1, with the geometry of each of the water molecule fixed into a rigid plane with C2v symmetry, having an OH distance of 0.9572 Å and a collapsed tetrahedral angle of 104.5. The potential between different water molecules i and j is given by 8 2 !12 !6 39 rc

ð5Þ

uc ðr1 , r2 , :::, rn Þ ¼

where the center of mass of the cluster (for equivalent masses), R, is given by R ¼

1 n rj n j¼1

∑

ð6Þ

The parameter rc is called the constraining radius with an accompanying constraining volume, most conveniently taken to be the volume of a sphere with vc = (4/3)πrc3, emphasizing that its center corresponds to the center of mass of the cluster. Because the introduction of the constraining potential is arbitrary, we take advantage of the inherent ambiguity to define the constraining potential with the center of mass of the oxygen atoms alone. As discussed elsewhere,31 the constraining potential imposes an internal pressure on the clusters. The choice of constraining radius is consequently of concern. It is important to choose rc to be sufficiently large that the properties of the clusters, at least at the lower temperatures, are not affected by the constraining potential. At the same time, if rc is chosen to be too large, the resulting simulation can have ergodicity problems.34 The specific values of the constraining radii used in this work are provided in section 3.1.2. 2.3. Cluster Partition Function. In addition to confining the clusters with a constraining potential, we ignore clustercluster interactions. Consequently, we can think about each cluster as an ideal-gas molecule, and to each cluster we can associate a molecular partition function. Classically, the molecular partition function for a cluster consisting of n indistinguishable molecules is defined by " !# Z ¥ n p2 1 j dp :::dpn exp β qn ðV , TÞ ¼ n!h3n ¥ 1 j ¼ 1 2m Z dr1 :::drn exp½βun ðr1 , :::, rn Þ ð7Þ

∑

V

where h is Planck’s constant, V is the macroscopic system volume, pi is the collective momentum of water molecule i in the cluster, ri represents the collective spacial coordinates of molecule i, and un is the interaction potential for the cluster. For an ideal gas consisting of clusters of a fixed size, the full, canonical partition function Q(T,V,N) can be written Q ðT, V , NÞ ¼

qn ðT, V ÞN N!

ð8Þ

where N represents the total number of clusters in the macroscopic system. Thermodynamic properties can be derived in the 4727

dx.doi.org/10.1021/jp201082p |J. Phys. Chem. B 2011, 115, 4725–4744

The Journal of Physical Chemistry B

ARTICLE

usual fashion from the partition function. Owing to the high dimensionality of the integrals required to evaluate the partition function, we rely on canonical averages using Monte Carlo methods to determine the thermodynamic properties of interest. 2.4. Gibbs Free Energy. A goal of this work is to determine the free energy changes associated with water cluster growth ðH2 OÞn þ H2 O f ðH2 OÞn þ 1

ð9Þ

and the addition of a hydrogen molecule to a water cluster ðH2 OÞn þ H2 f ðH2 OÞn H2

ð10Þ

2.4.1. Working Expressions. Because we are treating the clusters as ideal-gas particles, we begin with the standard-state ideal-gas expression for the chemical potential qn ðT, V Þ kB T μn ¼ RT ln ð11Þ V P where kB is the Boltzmann constant, the superscript denotes the standard state of 1 bar pressure, and R is the gas constant. Since the cluster is defined relative to its center of mass, it is convenient to transform the integral in configuration space to center-of-mass coordinates in eq 7 using the transformation n Y

drn ¼ n dR 3

j¼1

nY 1 j¼1

0

drj

ð12Þ

where the center of mass, R, and the relative coordinates, r0j are related by 0

rj ¼ rj R

j ¼ 1, 2, ... , n 1

ð13Þ

mole can be obtained from eq 17 and the expression ΔGn ¼ μn þ 1 μn μ1 "

! # P ðn þ 1Þ2 zn þ 1 ¼ RT ln kB T n3 zn

vn

where γ, the inverse thermal de Broglie wavelength cubed, is given by 2πmkB T 3=2 γ¼ ð15Þ h2 In eq 14 the integration over R yields the macroscopic volume, V, and the remaining integral is over the volume of the cluster with an interaction potential un relative to the clusters center of mass. Defining the configurational integral zn by Z Z 0 0 0 0 ::: exp½βun ðr1 , :::, rn 1 Þdr1 :::drn 1 ð16Þ zn ¼ vn

we obtain qn ðT, V Þ ¼ V

2πnm βh2

3=2

γn 1 n3=2 zn n!

ð17Þ

For the reaction given in eq 9 the Gibbs free energy change per

ð19Þ

Nearly identical expressions apply to the determination of the Gibbs free energy change for the reaction expressed in eq 10. To determine the change in the Gibbs free energy, we require a Monte Carlo approach for the calculation of the ratio of configurational integrals appearing in eq 19. 2.4.2. Thermodynamic Integration. Although numerous alternate techniques for calculating free energy changes exist, like perturbation methods23,41 and the nonequilibrium work methods of Jarzynski,42,43 for all of the applications considered here, we use a scaled Hamiltonian to carry out thermodynamic integration44 as developed by Kirkwood.23 In the Kirkwood formalism applied to water cluster growth, a dual-topological space is considered in which coexisting Hamiltonians of the initial and final states are connected through a so-called “alchemical transformation”45,46 via a coupling parameter λ. To evaluate the ratio of configurational integrals in eq 19, a simple choice for the dependence of the Hamiltonian on λ is HðλÞ ¼ ð1 λÞH0 þ λH1

ð20Þ

When λ = 0, H(λ) defines an initial system state and when λ = 1, H(λ) defines a final state. For the process expressed in eq 9 the initial state is a water cluster of size n and the final state is a water cluster of size n þ 1. To use the scaled Hamiltonian, we define a single-cluster analogue of the Helmholtz free energy, A, using the expression A ¼ kB T ln qn ðT, V Þ

3

The preproduct factor n in eq 12 is the Jacobian for the transformation between the primed center-of-mass coordinates and the unprimed coordinates. After integration of the momenta, eq 7 becomes Z γn 3 qn ¼ n dR n! V Z Z 0 0 0 0 ::: exp½βun ðr1 , :::, rn 1 Þ dr1 :::drn 1 ð14Þ

ð18Þ

ð21Þ

With the λ parametrization of the Hamiltonian, A becomes a function of λ, because AðλÞ ¼ kB T ln qn ðλÞ

ð22Þ

with the molecular partition function qn(λ) given by " !# Z ¥ n p2 1 j dp :::dpn exp β qn ðλÞ ¼ n!h3n ¥ 1 j ¼ 1 2m Z dr1 :::drn exp½ βun ðλ, r1 , :::, rn Þ

∑

ð23Þ

V

where the λ dependence of the Hamiltonian is apparent in the potential, un(λ), only. The temperature and volume dependence of qn(λ) have been ignored for notational convenience. Because the kinetic part of the Hamiltonian is independent of λ, and its contribution to the molecular partition function is analytic, the λ parametrization of the Hamiltonian in eq 20 can be replaced with the potential, un(λ), to give un ðλÞ ¼ ð1 λÞu0 þ λu1

ð24Þ

It is not difficult to verify that the derivative of A(λ) with respect to λ can be written as an ensemble average DAðλÞ Dun ðλÞ DAz ðλÞ ¼ ¼ ð25Þ Dλ V , T Dλ Dλ λ V, T where Az represents the configurational part of A. Integrating the 4728


The Journal of Physical Chemistry B lambda derivative of Az with respect to λ yields Z 1 Dun ðλ0 Þ ΔAz ¼ dλ0 0 0 Dλ 0 λ zn þ 1 ¼ RT ln zn

ARTICLE

ð26Þ

ð27Þ

where ΔAz = Az(λ=1) Az(λ=0). Equation 26 can be used directly to evaluate eq 19. 2.4.3. Singularities in the λ Integration. It is well-known47 that singularities in the Lennard-Jones potential at small distances lead to singularities in the integrand of eq 26. From scaled particle theory,48 one can show that in the case of the LennardJones interaction, the free energy derivative scales as λ3/4 as λ f 0, and this singularity is integrable, with a free energy that scales as λ1/4.49,50 However, the integrand of eq 26, and its subsequent integration, must be determined numerically in order to calculate the total free energy change. Consequently, although the singularity is integrable, the integrand must remain finite for numerical reasons. We have overcome the numerical difficulty associated with the singular behavior by modifying the parametrization of the Lennard-Jones potential energy function using the specific form introduced by Beutler et al.47

LJ

8 > > > > >
> > > > =

1 1 6 #2 > rij 6 > > > 2 rij > > 2 ð1 λÞ R > > LJ > ; : RLJ ð1 λÞ σ OO σ OO >

uij ¼ 4εOO λ "

ð28Þ In the denominator, RLJ is a positive constant and σOO and εOO are the Lennard-Jones parameters, shown in Table 1. It has been demonstrated47 that both of these initial states associated with the different values of RLJ (RLJ = 0 in the unmodified case) are equivalent in limλf0, rendering the free energy change independent of this choice. As discussed elsewhere,47 care needs to be taken in the choice of the parameter RLJ. We have found that by choosing RLJ to be in the range 0.01 e RLJ < 0.4, we are successfully able to evaluate numerically the integrals with respect to λ. Within this range of RLJ, no convergence problems in the derivatives of the potential have been encountered at small λ, and choosing RLJ = 0.05 has been found to be adequate for all of the water clusters studied. The singularity in the interaction between the charged sites of the Coloumbic part of the potential of eq 1 also gives rise to singularities in the integrand of the free energy. The Coulomb interactions can also be parametrized in a manner discussed by Beutler et al.47 qi qj uCij ¼ λ ð29Þ ½RC ð1 λÞ2 þ rij 2 1=2 in atomic units and RC is a positive constant. The Lennard-Jones repulsive term is expected to dominate the interactions at small λ, because the curvature of the Lennard-Jones potential energy function is much greater than that of the Coulombic function owing to the larger power of rij in denominator of the potential. With the dominance of the Lennard-Jones term at small λ, the inclusion of the Coulombic parametrization might seem overly cautious. However, in TIP4P water, the Coulombic interaction site is offset from the Lennard-Jones interaction site. Consequently, at

small λ, the Coulombic force may grow large at this interaction site before the Lennard-Jones repulsive term dominates at a different interaction site. As long as the value of RLJ is within an acceptable range, we have verified that setting RC = 1 is adequate for all of the water clusters studied. 2.4.4. Temperature Integration. For any physical or chemical process, the constant volume free energy change at any temperature, T, can be computed, in principle, by using the GibbsHelmholtz equation: D ΔA ΔE ¼ 2 ð30Þ DT T T V In eq 30, ΔE is the change in the thermodynamic internal energy function per cluster that includes both potential and kinetic energy contributions. To compute the free energy change for water cluster growth, the ensemble average of the derivative in eq 26 can be found using MMC with parallel tempering, at each temperature, for all λ. It is computationally more efficient, however, to do the λ integration at only one defined reference temperature, Tr, followed by a temperature integration. For the process given in eq 9, integration of eq 30 yields for the configurational part of the Helmholtz free energy change, Z T ΔAz ðTÞ ΔAz ðTr Þ ΔÆuæ 0 ¼ dT ð31Þ 02 T Tr Tr T To determine ΔAz(T), the free energy change at the temperature of interest, both terms on the right-hand side of eq 31 need to be calculated. The first term is the free energy change at some reference temperature, Tr, and is determined from eq 26. This λ integration is done numerically using a GaussLegendre quadrature scheme. Because the λ integration is performed at a single temperature, parallel tempering is inoperative. The reference temperature, Tr, is chosen high enough so that, in the absence of parallel tempering, the simulation is unlikely to suffer from any quasiergodicity problems. Very high temperatures can be in the dissociative region of the cluster. The dissociative region of temperature appears to offer an appealing choice for the reference temperature, since as the temperature increases, ΔAz approaches a noninteracting limit. Using such a limit, the λ integration becomes unnecessary, simplifying the calculation of ΔAz(T). Although such an approach has been used in previous work,30,51 we have found it to be most accurate to choose a reference temperature Tr, within a range of temperatures that are low enough to be considered in the “liquidlike” region of the cluster but high enough to sample sufficient phase space to keep the simulation ergodic. The second term on the right-hand side of eq 31 is the definite integral of the average energy difference between the different sized clusters divided by the square of the temperature with limits defined over the temperature range of interest. A plot of ΔÆuæ/T2 against T, shown in Figure 1, qualitatively illustrates the structure of the integrand for the reaction (H2O)4 þ H2O f (H2O)5. The upper curve marked by triangles shows the integrand diverging to more negative values as T f 0. At low temperatures, however, it is recognized that ΔÆuæ approaches its constant equipartition value, which implicates the T2 denominator of the integrand for this unwanted low-temperature behavior. Figure 1 also shows the integrand multiplied by T2, which transforms the functional behavior of the upper triangle curve into a curve, which is quite flat throughout the entire temperature range. This picture of the integrand with and without its denominator is suggestive of a 4729



Figure 1. Integrand of the temperature integration in eq 31 times a scaling factor, γ, plotted against T for the reaction (H2O)4 þ H2O f (H2O)5. The upper curve with data represented by triangles is the integrand with γ = 1 and diverges to more negative values at low temperatures. The lower curve with data represented by circles is the integrand with γ = T2 and is a smooth monotonic increasing function of the temperature.

Gaussian quadrature rule that can be adopted in a numerical integration scheme for this term. In Gaussian quadrature, the general treatment of numerical integration is changed from a single R b function integrand to an integration problem of the form af(x) w(x) dx. The weight function, w(x), can be introduced to handle singularities in the integrand. The flat behavior of the integrand on removal of the denominator in the temperature integration is suggestive of a generalized Gaussian quadrature scheme with 1/T2 as the weight function.52 By transferring this low temperature anomalous behavior to the weight function, the functional dependence of the integrand on the denominator has been effectively removed, resulting in the smooth, monotonic increasing function shown in Figure 1. A relatively flat curve of this type is very undemanding of a Gaussian quadrature, and the numerical method provided here is expected to require only a few Gaussian points and weights to yield accurate values for the integral in eq 31. To test the claim discussed in the previous paragraph, a reference temperature of Tr = 100 K has been chosen for the λ integration followed by a temperature integration using the proposed generalized Gaussian quadrature scheme for the 4-mer to 5-mer cluster reaction. A 12-point GaussLegendre quadrature is used in the λ integration at 100 K, using 2 107 MMC cycles, followed by temperature integrations with Tr = 100 and T equal to 0.5, 5, 10, 20, 30, 40, 50, and 80 K, using the same number of MMC cycles as the λ integration in the determination of ΔÆuæ. The points and weights used for the temperature integration have been generated from an Algol program53 that has been translated into Fortran. For the nine temperatures chosen to evaluate the free energy, the total number of evaluations of the average energy for the two-point quadrature is 2 8 = 16 (exclusive of the reference temperature, Tr, which is determined in the λ integration) plus an additional 11 temperatures for parallel tempering giving a total of 28 temperatures. As the number of quadrature points increases, the number of temperatures requiring energy evaluations also increases so that 43 temperatures have been used in the 4-point quadrature and 54 temperatures have been used in the 6-point quadrature. The different-point quadratures yield values of ΔAz(T) that are

ARTICLE

indistinguishable to within the statistical noise, confirming that a two-point quadrature over this range of temperatures is sufficient for convergence of the temperature integral. In further support of the proposed methods, a λ integration over a full range of temperatures from 1 to 350 K yields values of ΔAz that agree at each of the temperatures to at least 3 significant figures. In addition, it is observed that, in the harmonic limit, the value of ΔAz matches the expected value of ΔÆuæ, and this agreement lends credence to both the choice of reference temperature and the generalized quadrature rule proposed here; both are needed to determine the low-temperature value of ΔAz. 2.5. Steepest Descent Quench. The connection between the structure of the potential energy surface and the observed thermodynamic properties is an important goal of the current work. To facilitate an understanding of these connections, during the Monte Carlo simulations we often perform an inherent structure analysis54 to determine the isomer distribution associated with a given temperature. In the current context, an isomer is represented by a local (or global) minimum on the potential energy surface. To that end, we use conjugate gradient techniques to determine the nearest inherent structure from a given input configuration generated during the Monte Carlo walk. In the TIP4P model, the water molecules are treated as rigid and, traditionally, rigid-body rotations are accomplished by using rotation matrices constructed from an appropriate choice of an Euler angle rotation sequence. Unfortunately, the use of Euler angles introduces the problem of “gimbal lock”, or loss of one degree of rotational freedom owing to the degeneracy of the angles. It is the presence of this singularity in the use of Euler angles to rigid-body rotations that largely motivates the use of quaternions. Quaternions5557 replace the Euler sequence of concatenated rotations by a single rotation, which maps the reference and body frames in exactly the same way, but with a single rotation axis. The configuration of each monomer of a water cluster is expressed as a 1 6 dimensional p-vector for use in a conjugate gradient58 scheme. The p-vector of each monomer encodes its center of mass coordinates and the vector part, qB, of the quaternion q. The gradient of the potential, which is a function of the p-vector in quaternion space, is determined numerically using a quadratic approximation to the derivatives.59 An inherent difficulty in the use of numerical derivatives in the conjugate gradient scheme is the inaccuracy of the discovered minimum. Consequently, in a further effort to refine the minimum, the conjugate gradient output has been passed into a Monte Carlo program accepting only downhill values of the potential energy and quenched to absolute zero. After the conjugate gradient minimization has been enhanced with a Monte Carlo quench, a look at the Hessian is necessary to confirm that an inherent structure has indeed been found.

3. RESULTS 3.1. Pure Water Clusters. 3.1.1. Structures. The structures of the lowest energy TIP4P minima for (H2O)n with 2 e n e 20 can be found at the Cambridge Cluster Database.60 The simplest of these structures is the dimer, (H2O)2, a chain of two water monomers, and much experimental1 and theoretical6163 effort has been directed toward understanding the nature of the water dimer. Since hydrogen bonds are the dominant interaction between water molecules, the single hydrogen bond of the water dimer makes it the archetype for studying the hydrogen-bonded 4730



ARTICLE

Figure 2. Heat capacity curves of the dimer (lower left) and trimer (upper left) and tetramer (lower right) and pentamer (upper right).

networks in more complex water clusters. The global minimum structures for 3 e n e 5 are rings, with cluster sizes n g 6 having cages structures. 3.1.2. Heat Capacity. The potential energy contributions to the classical constant volume heat capacity per cluster, CV, for the various n-atom water clusters studied are calculated from the energy fluctuations given by the statistical expression CV ¼

Æu2 æ Æuæ2 kB T 2

ð32Þ

The kinetic energy makes a temperature-independent contribution to the heat capacity and is ignored in the current discussion. The heat capacity curves for n = 25 are shown in Figure 2. The displayed fluctuations represent two standard deviations, and if the error bars are absent, then it can be assumed their magnitude is less than the size of the iconic symbols used to distinguish different curves. The constraining radius, rc, from eq 5, has been chosen to be respectively rc = 4.0, 4.0, 5.0, 5.0, 6.5, 6.5, 8.0, 9.0, 10.0, 11.0, 12.0, 12.5, 13.0, 13.5, 14.5, 15.5, 16.5, 17.5, and 18.5 bohr for n = 220. For larger cluster sizes, rc increases so that the phase space volumes, albeit restricted, are large enough to allow freedom of rearrangement. As discussed previously, in all the Monte Carlo simulations performed, the constraining potential has been conveniently chosen to act on the oxygen atoms of each monomer alone. All heat capacity curves for n = 25 show libration-to-free-rotation (LTFR) peaks at low temperatures, and at temperatures beyond these lowtemperature peaks all four curves appear rather featureless. These LTFR peaks result from stereoisomers that can be interconverted by rotations. They are distinct from the melting peaks seen in larger clusters associated with phase changes that result in structural isomers with different hydrogen-bond connectivity patterns. To investigate the nature of the heat capacity curves of the water trimer, quench studies have been performed on configurations obtained from Metropolis Monte Carlo (MMC) simulations with parallel tempering using a quaternion modified conjugate gradient scheme followed by a downhill Monte Carlo scheme (see section 2.5). Configurations have been captured every 500 Monte Carlo passes on equilibrated structures for 2 105 MC cycles at each temperature, and then quenched to the nearest inherent structures. For the trimer, at temperatures of 150 K and below, shown as the top histogram in Figure 3, all configurations quench to the global minimum, with an energy of 69.99 kJ mol1. On all histograms shown in this work, U is

Figure 3. Histograms for the trimer isomer distribution at 150 K (top) and 300 K (bottom). The variable % p is the percent of the isomers of relative energy U compared to the total number of isomers analyzed in the inherent structure study.

the energy gap between higher energy inherent structures relative to the lowest energy isomer; i.e., the lowest energy isomer has been designated U = 0. All other isomers found at 150 K are sufficiently energetically inaccessible that we find no evidence for them on the scale of the number of Monte Carlo points included. At 300 K, 96% of the basins drain to the lowest energy isomer shown as the bottom histogram in Figure 3. The high-energy ring isomer of the trimer found at this temperature is nearly 18 kJ mol1 higher in energy than the lowest energy isomer. As illustrated in Figure 3, this ring isomer has the highest probability of occurrence among the high energy isomers. The configuration of the lowest energy isomer of the trimer along with its high-energy ring are illustrated in the left-hand and right-hand upper portion of Figure 4, respectively. One other local minimum has been found and is the open chain depicted as the bottom structure in Figure 4. The chain structure is 22.5 kJ mol1 higher in energy than the lowest energy isomer. As justified below, the low-temperature anomalies in the heat capacity curves of the dimer and trimer, with peaks at approximately 30 and 1.3 K, respectively, are a result of the onset of free rotations of the water molecules. Similar anomalies are observed in the heat capacity curves of the tetramer at approximately 40 K and the pentamer at 0.65 K. The low-temperature anomalies in the heat capacity curves for the trimer and pentamer are shown in the left- and right-hand sides of Figure 5, respectively. 4731



ARTICLE

In the first simulation, only rotational moves are allowed excluding translations, and the second simulation excludes rotations allowing only translational moves. In Figure 6 the upper curve is the heat capacity curve of dimer with its full complement of moves. Also displayed in Figure 6 are two curves generated from the limited moves discussed above. The middle curve is the heat capacity of the dimer constrained to rotations alone with a maximum that parallels the upper curve. It is observed that the increase in the rotational contribution coincides with the upper heat capacity curve of the dimer in the absence of translational moves. The lower heat capacity curve is limited to translational moves, and owing to anharmonic motion, is a monotonic increasing function of the temperature making little qualitative contribution to the heat capacity of the dimer in the temperature interval defined by a LTFR phase change. The appearance of a temperature region that marks the onset of free rotations has been observed in the heat capacity curves of (NH4Cl)n clusters,64 albeit at much higher temperatures than the peaks observed in the small water clusters presented here. The type of matching depicted in Figure 6 for the dimer has also been observed in the heat capacity curves of all the water ring systems. The water trimer is a more rigid structure than the water dimer bound by three nonlinear strained hydrogen bonds. As noted elsewhere, the alternation of free hydrogen atoms above and below the plane of the oxygen atoms, all with differing torsional angles, render this structure chiral,6567 as are all the oddmembered water rings. Figures 3.4 and 4.11 in Wales' book68

At the lowest temperatures considered in this work, the energy fluctuations are small, and the trimer is confined to its lowest energy configurations near the equipartition limit. As the temperature increases, the trimer begins to acquire sufficient energy to overcome the rotational barriers in configuration space separating the wells associated with different librations. As some molecules can occasionally extend their random walk to other wells from which they have been previously excluded, the energy fluctuations increase in order to surmount these librational barriers, and an associated increase in the heat capacity curve is observed. As the temperature increases, the rotations become free, the energy fluctuations decrease, and the heat capacity drops. The temperature at which this LTFR phase change occurs depends on both the magnitude of the rotational barriers between the wells and the number of degenerate wells that are accessible during this transition. These dependencies define energetic and entropic barriers, respectively, and each plays a key role in determining temperature-dependent structural changes. In qualitative support of this rotational transition from hindered to unhindered motions, two further simulations have been performed in which the Monte Carlo move strategy has been decomposed into separate translational and rotational parts.

Figure 4. The three discovered trimer isomers at 300 K. The global minimum is the upper left-hand structure, the upper right-hand structure is the highest energy ring, and the lower structure is the chain. All structures displayed in this work have been generated with xmakemol (http://www.nongnu.org/xmakemol), with oxygen atoms represented by red spheres, hydrogen atoms represented by white spheres, and dashed lines representing the hydrogen bonds.

Figure 6. (Upper curve) Heat capacity curve for the dimer for temperatures below 100 K. (Middle curve) Heat capacity curve for the dimer computed in the absence of translations. (Lower curve) Heat capacity computed with rotations excluded.

Figure 5. Low-temperature heat capacity curves of the trimer and pentamer. 4732



ARTICLE

Figure 7. Constant volume heat capacity curves for (H2O)n with n = 6 left curve and n = 11 right curve.

illustrates the energy differences associated with the degenerate flip rearrangements of the water trimer. The position of two adjacent hydrogen atoms on the same side of the ring makes the water trimer a frustrated structure, which results in easy unhindered librational motions with lower rotational energy barriers than the dimer.69,70 The heat capacity curves shown in Figure 2 for the trimer and pentamer are strikingly similar, as are the heat capacity curves for the dimer and tetramer. These similarities can be attributed to the chiral nature of the odd-membered rings and the achiral behavior of the even-membered rings. The difference between (H2O)n for n = 2,4 and n = 3,5 is the geometrical frustration of the chiral odd-membered rings giving a lower temperature anomaly than the hindered (unfrustrated) rotations of the achiral structures. There is an apparent nonadditivity of the translational and rotational contributions to the heat capacity; i.e., the sum of the computed translational and rotational heat capacities does not give the total heat capacity of the system. This nonadditive behavior contrasts with the observation in ammonium chloride64 and arises artificially, because for water the method of making rotational moves alone includes some translational contributions.71 Such rotational/vibrational coupling for rotational moves in ammonium chloride are absent owing to the spherical symmetry of the gas-phase ammonium ion (lacking in water). With the exception of n = 6 and n = 11, the constant volume heat capacity curves as a function of temperature for water clusters for 6 e n e 20 are similar to previously published heat capacities of (H2O)8.8,7275 All the clusters in this size range show a melting peak significantly lower than their bulk counterparts (the melting point of TIP4P bulk water reported by Vega and Abascal76 is 232 K), with the heat capacity rising from its equipartition value at T = 0 until the temperature approaches the cluster dissociative region. The curves for n = 6 and n = 11, displayed in Figure 7, illustrate some of the unique features of these exceptions, which are absent from the other sized water clusters in this range. For other sized clusters, detailed curves and discussions can be found elsewhere.71 Additionally, in section 3.2.1 in the discussion of the heat capacities of impure water clusters, the heat capacity curves for (H2O)10 through (H2O)20 are displayed as dashed lines in Figures 2022. The hexamer shows a low-temperature LFTR peak at 5 K in its heat capacity curve as well a melting peak associated with the onset of structural isomerization at 93 K. Other larger sized clusters display only a melting peak. The 11-mer, (H2O)11, has a shoulder in its heat capacity curve around 91 K, and the curve is

Figure 8. Histograms for the isomer distributions of (H2O)11 at 1 , 9 , 91, and 207 K from top to bottom, respectively.

extremely broad beyond its peak temperature at 225 K. To investigate the origin of the anomalies in the heat capacity curve of the 11-mer, quench studies of the configurations obtained from simulations at temperatures of 1, 9, 90, and 207 K have been performed. The histograms in Figure 8 show the isomer distributions for the 11-mer at these four temperatures. At each temperature, only isomers with energies less than 20 kJ mol1 above the lowest energy isomer are shown. The top histogram illustrates that at 1 K all configurations quench to the global minimum. In this case, the lowest energy isomer at U = 0 corresponds to an energy of 431.5 kJ mol1. At 9 K, the second histogram from the top in Figure 8, a nearly isoenergic pair of isomers have been found with the high-energy isomer only 0.17 kJ mol1 higher in energy than the lowest energy isomer. The structures for this near isoenergic pair of isomers for the 11mer, illustrated in Figure 9, show different hydrogen-bond bridging patterns between the water monomer at the apex of the pentagonal ring with the monomer of the recessed tetrameric ring. Differing hydrogen-bond bridges appear to induce slightly different dipole orientations (see section 3.1.4), which could explain the energy gap in this pair. Many of the observed isomers for the 11-mer show these alternating bridge patterns with differing degrees of distortion between the rings. At 91 K, the third histogram from the top in Figure 8, numerous isomerizations in this “premelting” region of the 11-mer have been identified. These isomerizations include bridged cubes and 4733



Figure 9. Nearly isoenergetic pair of isomers of (H2O)11 discovered by a steepest descent quench at 9 K.

ARTICLE

followed by the 4-point generalized Gaussian quadrature scheme for the temperature integration discussed in section 2.4.4. The standard-state Gibbs free energy change is then determined from eq 19. For each cluster, ΔAz has been determined for a series of eight temperatures, 0.5, 5, 10, 20, 30, 40, 50, and 80 K, exclusive of the reference temperature (100 K). Additional temperatures have been added to the temperature points (and weights) used in the numerical evaluation of the integral in eq 31 to ensure an appropriate implementation of parallel tempering. For the simulation of the trimer and pentamer, a 5 107 MMC cycle equilibration is followed by 1 108 cycles for data accumulation. A 3 107 equilibration period followed by 6 107 cycles of accumulation have been performed for dimer and tetramer. For each of these simulations, parallel tempering is used with a temperature grid having a probability of acceptance of at least 30% for each exchange. Using our initial assumption that the water clusters behave ideally, the temperature-dependent enthalpy change, ΔHn(T), is for the related to the molar internal energy change, ΔE(T), n reactions given by eq 9, by ΔHn ðTÞ ¼ ΔEn ðTÞ RT

ð33Þ

The potential energy change is determined from the caloric curves computed during the Monte Carlo simulations of each water cluster. The entropy change, ΔS, for each process is obtained from the thermodynamic expression ΔSn ¼ ðΔHn ΔGn Þ=T Figure 10. Isomers of (H2O)11 found by steepest descent quench in the “premelting” region of the temperature.

bridged pentagonal prisms. These types of structural changes are illustrated in Figure 10. The configuration in the left-hand side of Figure 10 displays a pentameric face attached to the edge of a cube, which is quite different from the type of bridging seen in Figure 9. The right-hand side of Figure 10 shows bridging hydrogen bonds across a pentagonal prism. The distinctly different configurations of these isomers are relatively close in energy, differing by no more than 1.0 kJ mol1. These isomers are approximately 3 kJ mol1 higher in energy than the lowest energy isomer. The low-temperature anomaly in the heat capacity curve of the 11-mer is attributed to frequent visits between a large number of nearly isoenergic isomers with different structures. The lack of similarity between these topologies is sufficiently distinct to generate energy fluctuations that cause the observed anomaly in this temperature region of the heat capacity curve. The broadening of the peak in the melting region of (H2O)11 is attributed to the presence of more isomeric forms. The isomer spread in this region of the 11-mer is shown in the bottom histogram of Figure 8. 3.1.3. Gibbs Free Energy. The standard Gibbs free energy changes, ΔGn, for the reactions expressed in eq 9 with n = 2 e n e 5 are shown in Figure 11. To determine the free energy change, ΔG, n for the dimer and water rings a reference temperature of Tr = 100 K has been chosen for both the λ and temperature integration, which is a temperature well below the dissociative region of these small clusters, but a temperature that is sufficiently high to prevent entrapment in local basins on the potential energy surface. The free energy change for each of these clusters is determined by first calculating ΔAz(Tr) numerically using a 12-point GaussLegendre quadrature for the λ integration

ð34Þ

Each of the panels in Figure 11 shows the temperature dependencies for each of these thermodynamic state functions over the indicated temperature range. The general trends in ΔG, n ΔH, n and TΔSn for each of these processes are very similar for each of these small water cluster systems, and these trends are observed to persist throughout the entire series of processes for the cage clusters as well.71 The limiting behavior of ΔGn at both high and low temperatures can be understood qualitatively by the general expression for the standard free energy change decomposed into its enthalpic and entropic parts. ΔGn ¼ ΔHn TΔSn

ð35Þ

In eq 35, as T f 0, ΔGn f ΔHn, a limit that is clearly illustrated in Figure 11 for the dimer and water ring structures. Since the number of reactant particles in the reaction mixture is decreasing, the entropy change for the reaction of water cluster growth is negative. The free energy changes for the reactions given in eq 9 for the cage clusters with 6 e n e 20 are qualitatively identical71 to those we have already examined for n e 2 e 5. Consequently, rather than exploring the temperature dependence of these larger clusters, we now examine the standard free energies as a function of cluster size for several fixed temperatures. For water clusters sizes beyond the water rings, higher reference temperatures have been chosen to bracket temperature regions that include features of the cluster transitions that are physically significant. For the hexamer, the temperature integrations for ΔAz have been determined with a reference temperature of Tr = 220 K at a series of 12 other temperatures, and for all other water cages up to n = 20, a reference temperature of Tr = 320 has been chosen. Additionally, ΔAz has been determined at series of 17 other temperatures. Again, the free energy change for each of these cage clusters is determined by first calculating ΔAz(Tr) 4734



ARTICLE

Figure 11. In all the panels, ΔX, expressed in kJ mol1, is either ΔG, n the middle (circles) curve, or ΔH, n the lower (up triangles) curve, or TΔS, n the upper (diamonds) curve. For each of the processes indicated, ΔHn is nearly constant over the 0.5100 K temperature range. The upper-left panel is for the n = 23 process, the upper-right panel represents n = 34, and the lower panel represents n = 45.

numerically using a 12-point GaussLegendre quadrature for the λ integration followed by a 4-point generalized Gaussian quadrature scheme for the temperature integration. For the simulation of these water cages, a 4 107 MMC cycle equilibration is followed by the same number of cycles for data accumulation. A graph of ΔGn as a function of n is plotted in Figure 12. Two outstanding local minima occur for n = 7 and n = 11 associated with the processes (H2O)7 þ H2O f (H2O)8 and (H2O)11 þ H2O f (H2O)12, respectively. The structures of the 7-mer and 11-mer can be viewed as a cubiclike structure with a vacancy at a corner. Both of the reactions, n = 7 and n = 11, fill this vacancy with a water monomer to give a cube for n = 8 and a fused cube for n = 12. Figure 12 also displays the effect of raising the temperature of each of these processes. As the temperature increases from 5 K, the free energy change for each of these processes is progressively diminished as the system becomes more fluid, and at 220 K any vestiges of the defining structural characteristics of the systems have vanished. In Figures 13 and 14, the temperature behavior of the enthalpic and temperature-weighted entropic components (i.e., TΔS) n of the free energy respectively are shown as a function of n. While the enthalpy changes are nearly temperature independent, the temperature-weighted entropy changes are the antithesis of those observed for the free energy. The temperature-weighted entropy change alone, however, is no larger than 0.13 kJ mol1 over the entire cluster series considered here and is unimportant in determining relative stabilities. In an attempt to gain some insight concerning the origin of the relative stabilities of the clusters as measured by their free energy

Figure 12. Free energy change, ΔGn, for the processes (H2O)n þ H2O f (H2O)nþ1 at T = 5 K, T = 100 K, T = 160 K, and T = 220 K. Each process is identified by a point given by the value of n, and the points are connected by straight lines. As the temperature increases, the structure is progressively “washed-out” of each process.

difference as a function of cluster size, we define the notion of the number of hydrogen-bonded interactions in a given cluster structure. The hydrogen-bonded pairwise monomer interactions in this study is defined to have an interoxygen distance rOO e 5.5 bohr and θ g 120, where θ is OH 3 3 3 O angle. The OO distance between neighboring monomers and its strength is given by the pairwise TIP4P potential for the OH 3 3 3 O interaction. For TIP4P water clusters, at the putative global minima listed in the Cambrdige Cluster Database,60 the number of OO nearest neighbors is the same as the number of hydrogen bonds, which provides a convenient method of 4735



Figure 13. Enthalpy change, ΔH, n for the processes (H2O)n þ H2O f (H2O)nþ1 at T = 5 K, T = 100 K, T = 160 K, and T = 220 K. Each process is identified by a point given by the value of n, and the points are connected by straight lines. The enthalpy changes for all of these processes are observed to be virtually temperature independent and match the behavior of ΔGn at low temperatures.

ARTICLE

Figure 15. Energy per hydrogen bond versus temperature.

Table 3. Classification Scheme for Pairwise Monomer Interaction (PMI) Strength Relative to the Single Hydrogen Bond of the Dimer

Figure 14. Temperature-weighted entropy change, TΔSn, for the processes (H2O)n þ H2O f (H2O)nþ1 at T = 5 K, T = 100 K, T = 160 K, and T = 220 K. Each process is identified by a point given by the value of n, and the points are connected by straight lines.

hydrogen bond identification. All pairwise interactions that are OO nearest neighbors are hydrogen-bonded pairs, and all pairs of monomers that are not OO nearest neighbors are nonhydrogen-bonded pairs. Using the foregoing criterion for determining the number of hydrogen bonds in the putative global minimum structures for each cluster size, the first two columns of Table 2 list the number of hydrogen bonds for each cluster size. The energy per hydrogen bond, UHb, for all of the cage clusters, 6 e n e 20, is shown in Figure 15 as a function of temperature (over a limited temperature range for clarity). The average energy per hydrogen bond for the clusters, displayed in Figure 15, offers little insight into the fine structure observed in Figure 13 for these enthalpy changes (and Figure 12 for the lowtemperature free energy changes). To acquire the acumen necessary to understand the nature of these changes, a closer look at the collective orientations of the monomer dipoles provides an important component of a useful model. 3.1.4. Model for the Energetics of Cluster Growth. As discussed in section 3.1.3, it is not possible to rationalize the fine structure in the enthalpy change as a function of cluster size (or the free energy change as a function of cluster size at low temperature) by counting the number of hydrogen-bonded near-neighbor interactions. In other systems77 modeled by pair

strong

90% < PMI < 100%

moderate

70% < PMI < 90%

weak

30% < PMI < 70%

very weak

PMI < 30%

potentials, such correlations between the number of near neighbors and the enthalpy change have been useful, and we now explore those aspects of the energetics of water clusters that help us understand the trends observed in Figure 13. A principal difference between the observations for water clusters and other systems is the long-range nature of the Coulombic interactions, and we now show that we can model the energetics in a way that can help us understand the observed fine structure. In section 3.1.3 we have introduced a distance criterion for deciding whether two water molecules are directly hydrogen bonded to each other. Using that criterion, we can then associate a pair energy for each interaction in the cluster. We normalize each of the pair energies with respect to the dimer energy and the strength of the resulting pairwise monomer interactions (PMI) are then classified as either “strong”, “moderate”, “weak”, or “very weak” according to the scheme shown in Table 3. We find that the energies of the associated hydrogen bonds vary significantly, and this observed variation in hydrogenbonded PMI strengths, shown in Table 2, correlates poorly with the OO distances, which vary only slightly from approximately 5 bohr. An understanding of the rather large disparity in the hydrogen-bonded PMI strengths over near-equidistant oxygen centers can be obtained by examining the orientation of the monomer dipoles. For every cluster, each TIP4P monomer dipole has been programmed into a molecular visualization program (see the caption to Figure 4), and for each pair of monomers involved in hydrogen bonding the classified PMI strengths for the pair is identified. A direct correspondence between the orientation of the dipoles between each pair of monomers and the observed hydrogen-bonded PMI strength leads to the following two postulates. Postulate 1. If the dipole orientations on neighboring monomers are cis relative to the cylindrical symmetry along the internuclear axis of the oxygen atoms defining the hydrogen 4736



ARTICLE

Figure 16. Two views of the S4 global energy minimum of (H2O)8. Figure 17. Two views of the D2d local energy minimum of (H2O)8.

bond, the pairwise monomer interactions (PMI) are relatively weak. The more parallel the cis dipoles, the weaker the PMI. Postulate 2. If the dipole orientations on neighboring monomers are trans relative to the cylindrical symmetry along the internuclear axis of the oxygen atoms defining the hydrogen bond, the PMI are relatively strong. The more parallel the trans dipoles, the stronger the PMI. The last two columns in the Table 2 give the total number of trans and cis dipoles associated with hydrogen bonds in each water cluster. All hydrogen-bonded PMIs classified as either “strong” or “moderate” are trans-type PMI and any classified as “weak” or “very weak” are cis-type PMI. The dipole orientations on monomers involved in hydrogen bonding, however, change drastically with changes in cluster structure and strongly influence the energetic contributions to the total energy. If the dipole orientations on neighboring monomers connected by hydrogen bonds are cis-type, the pairwise interactions are considerably less attractive. Conversely, if the dipoles are trans-type, the pairwise interactions are much stronger. To illustrate the model, we now consider the various energy contributions to the octamer, (H2O)8. The global minimum energy structure of the octamer has S4 symmetry, and two views of the S4 structure are shown in Figure 16. The view on the lefthand side of Figure 16 illustrates the notion of homodromy.19 The four-membered rings in the front and recessed faces of the lefthand view have the hydrogen bonds in both rings connected in the same sense, and the homodromy is deemed parallel. Shortly, the parallel homodromy found in the S4 isomer is contrasted with the antiparallel homodromy of the isomer having D2d symmetry. In Figure 16 the arrows indicate the directions of the dipole moments of all the water molecules, and in this parallel view, there are four cis-type dipoles between monomers connecting these two rings and another four trans-type dipoles between vertical monomers on the same tetrameric face. The right-hand view of Figure 16 shows the orientation of the dipoles through the square face rotated by 90 from the left-hand view. The trans orientation of the remaining four hydrogen-bonded dipoles is apparent in this view. Only slightly higher in energy than the S4 global minimum energy structure of (H2O)8 is the structure with D2d symmetry depicted in Figure 17. The sense of the hydrogen bonds in the two rings depicted on the left-hand view in Figure 17 is antiparallel, and the structure has antiparallel homodromy. As in Figure 16, the right-hand view in Figure 17 is rotated from the left-hand view by 90. As is clear in Figure 17, the dipole interactions are all trans in the D2d isomer. A comparison between the hydrogen-bond strengths in the S4 and D2d isomers of (H2O)8 is given in Table 4. If the relative

Table 4. Distribution of Hydrogen-Bonded PMI Strengths for the S4 and D2d Isomers Shown in Figures 16 and 17 isomer n hydrogen bonds strong moderate weak very weak trans cis S4

8

12

8

0

4

0

8

4

D2d

8

12

12

0

0

0

12

0

energies of the two isomers were determined by the hydrogen bonds alone, we might expect the D2d isomer with all-trans orientations to be lower in energy than the S4 isomer containing both cis and trans dipole orientations. However, the S4 isomer is lower in energy, and the energy order implies the need to consider the contributions from the non-hydrogen-bonded interactions to understand the observation. The total number of pair interactions in a cluster of size n is n(n 1)/2. For the octamer the total number of pair interactions is 28. Using our definition of hydrogen-bonded interactions, both the S4 and D2d isomers have 12 pair interactions that arise from hydrogen bonding. To understand the energy ordering of the two isomers, we must consider both hydrogen-bonded pair energies and the contributions from the pair interactions that are not hydrogen bonded. By examining each of the pair interactions that are not hydrogen bonded, we find some are repulsive (0). In Table 5 we define Erep to be the sum of the repulsive non-hydrogen-bonded pair interactions, Eatt to be the sum of the attractive non-hydrogenbonded pair interactions, EH to be the sum of the hydrogenbonded PMI, and Etot to be the total energy. The data in Table 5 make clear that while the hydrogen-bonded energy in the D2d isomer is lower than the hydrogen-bonded energy for the S4 isomer, as expected from the dipole orientations, the repulsive energy in D2d is large compared to the near-zero repulsive energy in the S4 isomer. The difference in repulsive energy, Erep, between these isomers is largely a result of the repulsively coherent dipoles observed along the face diagonals of the D2d cube in Figure 17. Such strongly directed repulsions are absent in any view of the S4 isomer. The small energy difference between the two isomers results from a complicated mix of contributions. We can continue this analysis for other size clusters as well. For the putative global minimum energy structures, the difference in the total PMI energy for those monomers connected by hydrogen bonds as a function of cluster size is plotted in Figure 18 with the enthalpy curve of Figure 13 superimposed for comparison. For cluster processes with n = 2, 5, 7, 12, and 15, the difference in hydrogen-bonded PMI energy closely mirrors the enthalpy change; in other cases the correlations with the enthalpy changes 4737



ARTICLE

Table 5. Distribution of PMI Energy in the Non-HydrogenBonded and Hydrogen-Bonded Interactions for the D2d and S4 Isomers at 1 Ka cluster

a

Erep

Eatt

EH

Etot

D2d

40.47

48.89

296.9

305.3

S4

≈0.0

40.46

265.0

305.5

All energies are reported in kJ mol1.

Figure 19. Difference in PMI energies as a function of cluster size, n, for the processes (H2O)n þ H2O f (H2O)nþ1. The open circles identify near neighbor, hydrogen-bonded PMI energies (Figure 18), and the filled circles identify the non-near neighbor (non-hydrogen bonded) PMI processes for each value of n. The points are connected by straight lines. The enthalpy changes from Figure 13 are displayed (the dashed line) for comparison.

Figure 18. Difference in hydrogen-bonded PMI energies as a function of cluster size, n, for the processes (H2O)n þ H2O f (H2O)nþ1. Each process is identified by an open circle given by the value of n and the points are connected by straight lines. The dashed line is the enthalpy change taken from Figure 13 for comparison.

are poor, and the non-hydrogen-bonded interactions must explicate the observed energy difference. In the context of this heuristic, the “mixed interactions” energy, which includes both attractive and repulsive components, is defined as the total nonhydrogen-bonded PMI. This “non-near-neighbor” energy is equal to the change in “mixed interaction” energy, denoted , and can be expressed as ΔHnon-Hb ¼ ΔHn ΔHHb ΔHnon-Hb

ð36Þ

= Δ[cis-type where the change in hydrogen-bonded PMI, ΔHHb PMI þ trans-type PMI]. To account for this non-near-neighbor energy, all the non-hydrogen-bonded interactions between pairs of monomers have been determined for each cluster. The differences between these interactions for contiguous clusters along with the curves of Figure 18 are shown in Figure 19. The sum of the energies represented by the filled and open circles, for each value of n, gives the enthalpy change for the process. Figure 19 emphasizes the importance of these “mixed” interactions, defined as those pairwise interactions not associated with hydrogen bonding. These mixed interactions compose a large part of the observed enthalpy change for many of the clusters in this series. This analysis of the contributions to the energy changes can be applied to a number of cluster sizes in the series providing additional insights. The detailed discussion of the energetics of the clusters can be found elsewhere.71 3.2. Water Clusters with Molecular Hydrogen. 3.2.1. Heat Capacities. The constant volume heat capacity curves as a function of temperature for water clusters “doped” with hydrogen (H2O)nH2 ranging in size from 10 e n e 20 are shown in Figures 2022 as lines connecting data points that are

represented by circles. All of these curves have been determined over a temperature range of 1320 K in a manner similar to those for the pure water clusters. For each cluster, a 6 107 MC cycle equilibration period has been followed by the same number of MC cycles of data accumulation. The constraining radii used in the mixed clusters considered here are respectively 11.0, 12.0, 13.0, 14.0, 14.5, 15.0, 16.5, 17.0, 18.0, 19.0, and 19.0 bohr for n = 1020. This range of constraining radii is typically 1.01.5 bohr larger, for the same value of n, than the values used for the pure water clusters. The heat capacity curves for each of the underlying pure water clusters71 are superimposed for comparison (the dashed curves). In the main, all the curves agree, within statistical fluctuations, with the pure water clusters, except at the extremes of temperature. At low temperatures, all of the heat capacity curves for the impure water clusters show a low-temperature anomaly, which is absent in their pure cluster counterparts. In addition, these impure cluster heat capacity curves display more internal degrees of freedom in the equipartition limit. As a general rule, for classical systems the more internal degrees of freedom, the greater the heat capacity. In the harmonic limit, for a cluster with Nw water molecules and Nh molecules of hydrogen, [3(3Nw þ 2Nh) 6]/2 quadratic terms are introduced as potential energy into the Hamiltonian. For Nw water molecules and Nh hydrogen molecules, (3Nw þ Nh)/2 modes have been “frozen-out” from the assumption of rigidity for the nonlinear water molecules and linear hydrogen molecules. Consequently, for an impure water cluster with only one hydrogen molecule (i.e., Nh = 1), the heat capacity should approach the following equipartition limit lim

Tf0

CV 1 ¼ 3Nw 2 kB

ð37Þ

For example, for the heat capacity curve of the (H2O)10H2, shown in the upper left-hand corner of Figure 20, the lowtemperature heat capacity approaches the equipartition limit of 29.5kB. This value is 2.5kB higher than the associated lowtemperature heat capacity of (H2O)10, which has an equipartition value of 27kB. This disparity in heat capacities between pure 4738



ARTICLE

Figure 20. Constant volume heat capacity curves for (H2O)nH2 with n = 10 (top left), n = 11 (top right), n = 12 (bottom left), and n = 13 (bottom right). The dashed curve represents the heat capacity for the same sized pure water clusters.

Figure 21. Constant volume heat capacity curves for (H2O)nH2 with n = 14 (top left), n = 15 (top right), n = 16 (bottom left), and n = 17 (bottom right). The dashed curve represents the heat capacity for the same sized pure water clusters.

and impure clusters at the equipartition limit, which illustrates the higher degrees of freedom present in this heterogeneous

cluster as T f 0, can be observed in Figures 2022 for every cluster studied in this work. 4739



ARTICLE

Figure 22. Constant volume heat capacity curves for (H2O)nH2 with n = 18 (top left), n = 19 (top right), and n = 20 (bottom). The dashed curve represents the heat capacity for the same sized pure water clusters.

The equilibration of randomly initialized impure water clusters has been found to be very expensive computationally, especially for the larger sized clusters. Consequently, the initial structures in each simulation have been taken to be the coordinates of the global minimum for each of the water clusters and additional coordinates have been chosen that place the hydrogen molecule in close proximity, but “outside” of the water cage. The discovered putative global minimum configuration for (H2O)10H2 is shown in Figure 23, with the hydrogen molecule outside the main body of the water cluster. This equilibrium configuration is found to be of lower energy than any choice with the hydrogen molecule “inside” the main body of the water cluster. Further, for the inside initializations of these heterogeneous clusters, the hydrogen molecule never remains sequestered inside the water cage. For example, an inside random initialization of the cluster, (H2O)10H2, reveals a pentagonal prism with an outside hydrogen as the low-energy structure. In this case, and for other cases considered, the hydrogen remains outside the water skeleton throughout the entire range of temperatures considered. To investigate the low-temperature anomalies in these heat capacity curves (Figures 2022), additional simulations have been performed on (H2O)10H2, which has been chosen as a prototype for the low-temperature anomalies present in all the clusters considered in this work. For each of these additional simulations, various combinations of translational and rotational motions for the hydrogen solute and the water monomers have been performed. The heat capacity curve for the simulation of (H2O)10H2 in which the impurity is allowed to translate and rotate, but all comparable motions of the water monomers have been “frozen”, is shown as the lower curve in Figure 24. The upper curve is identical to the heat capacity curve shown in the

Figure 23. Putative global minimum structure of (H2O)10H2.

upper left-hand corner of Figure 20. The low-temperature anomaly in the heat capacity curve of this cluster at 10 K (upper curve) is seen to match the peak of the lower curve, providing qualitative evidence that the low-temperature heat capacity anomaly is a result of the hindered rotations and translations of the hydrogen molecule. It must be emphasized that both rotational and translational moves are required to reproduce the lower curve in Figure 24. The heat capacity anomaly for translational or rotational moves alone occurs at temperatures about 20 K higher than the anomaly when both motions are present. It appears that the presence of the additional mode lowers the energy barriers for either motion considered separately. The lower curve in Figure 24 rises to a maximum at about 10 K and then decreases. This rise is coincident with intermittent jumps to degenerate orientations and marks the 4740


The Journal of Physical Chemistry B onset of free rotations and translations. At this temperature, both of these motions are hindered by the presence of a relatively large hydrogen-bonded water skeleton. Beyond 10 K, the motions of the hydrogen molecule are more representative of free rotations and translations, rather than jumps, and with more accessible phase space, the energy fluctuations decrease. This decrease in energy fluctuations is accompanied by a decrease in the rotational and translational contributions of the hydrogen solute to the heat capacity. In the equipartition limit, for this restricted motion, CV/ kB = 2.5, consistent with the two rotational and three translational degrees of freedom for a linear molecule. At high temperatures, the

Figure 24. (Upper curve) Constant volume heat capacity curve for (H2O)10H2 as depicted in the upper, left-hand corner of Figure 20; (lower curve) heat capacity curve of this “mixed” cluster, in which the only moves that have been allowed are the translations and rotations of molecular hydrogen. All comparable motions of the water molecules are excluded when generating the data represented in the lower curve.

ARTICLE

impure clusters (Figures 2022) show a slightly higher heat capacity. Presumably, in the melting region, the presence of hydrogen results in a greater variety of isomers with more degrees of freedom than pure water clusters. 3.2.2. Gibbs Free Energy Change. As we have seen for pure water clusters, we find the qualitative behavior of ΔGn, ΔHn, and TΔSn as a function of temperature for the reaction expressed in eq 10 to be the same for all n. The changes for these thermodynamic quantities for 10 e n e 13 are shown in Figure 25. Because of the similarity of the curves, we do not show the data for n e 14 e 20.71 For the entire set of these impure water clusters, the temperature integrations for ΔAz have been determined at series of 17 temperatures, ranging from 1 to 300 K exclusive of a chosen reference temperature of Tr = 320 K. All of the curves shown in Figure 25 have been truncated at 50 K to render the low-temperature behavior visible. At higher temperatures there is no change in the general behavior. In a manner completely analogous to pure water clusters, the free energy changes for each of these clusters have been determined by first calculating ΔAz(Tr) numerically using a 12-point GaussLegendre quadrature for the λ integration followed by a 4-point generalized Gaussian quadrature scheme for the temperature integration as discussed in section 2.4.4. In this case, however, the Kirkwood coupling parameter, λ, targets the molecular hydrogen as the initial, homogeneous, n-monomer water cluster transforms into a heterogeneous (H2O)nH2 cluster through a series of intermediates. For each of these impure water clusters the form of the scaled Hamiltonian given in eq 20 is identical to that used for pure water clusters. For the simulation of these hydrogen “doped” water cages, a total of 2.0 108 MMC cycles have been run for equilibration and data accumulation.

Figure 25. For the addition of H2 to (H2O)10 (upper left), (H2O)11 (upper right), (H2O)12 (lower left), and (H2O)13 (lower right), ΔX = ΔG (circles), ΔH (up triangles), or TΔS (diamonds), expressed in kJ mol1. For each of these exothermic processes, ΔH is nearly constant over the entire temperature range. 4741



Figure 26. Free energy change, ΔGn for the processes (H2O)n þ H2 f (H2O)nH2 at T = 5 K, T = 100 K, T = 160 K, and T = 220 K. Each process is identified by a point given by the value of n, and the points are connected by straight lines.

The trends in ΔGn, ΔHn, and TΔSn are similar to those of the pure water clusters for the same reasons given in section 3.1.3. In particular, at the low-temperature limit where the equipartition theorem applies, ΔG = ΔH. The free energy change is about 3 kJ mol1 at low temperatures. This free energy difference is substantially smaller than those observed in pure water clusters, and the lack of variation with increasing cluster size, n, betokens minimal changes in the water component of the cluster as hydrogen is added. A graph of the standard free energy change with increasing cluster size, ΔG, n for each of the hydrogen cluster reactions in eq 10, is plotted in Figure 26 for a range of temperatures. It can be observed that ΔGn is nearly independent of n for all T. This behavior is in stark contrast to the free energy behavior of pure water clusters, shown in Figure 12. The free energy components, ΔHn and TΔSn, are plotted in Figures 27 and 28, respectively. The enthalpy change is seen to be independent of temperature with little structural change over the range of cluster sizes considered. Again, this enthalpy behavior contrasts with the enthalpy behavior observed in the case of water cluster growth, shown in Figure 13. The underlying reason for the seemingly soporific behavior of these impure water clusters can be understood by pondering the structure depicted in Figure 23. Figure 23 (with n = 10) is representative of the product species, (H2O)nH2, for all of the reactions in eq 10. The pattern occurs because the hydrogen molecule remains “outside” the water cage for every n, and the addition of hydrogen does not modify the initial structure of the water cluster without hydrogen. Consequently, the only interactions observed in these clusters are between the hydrogen molecules and the skeletal moieties of the water cage. As such, these interactions are virtually independent of cluster size, since they depend only on the availability of the dangling hydrogen atoms on the surface of the water cages. Although larger clusters offer a greater number of these outward-pointing surface hydrogen atoms, the hydrogen molecule is only interacting with one of these moieties and remains insensitive to their number.

4. SUMMARY AND DISCUSSION In this work, we have investigated the thermodynamic properties of water clusters modeled with a TIP4P potential with and without the addition of molecular hydrogen. We have focused on those properties that relate to the thermodynamic stability of the

ARTICLE

Figure 27. ΔHn for the processes (H2O)n þ H2 f (H2O)n(H2) at T = 5 K, T = 100 K, T = 160 K, and T = 220 K. Each process is identified by a point given by the value of n, and the points are connected by straight lines.

Figure 28. TΔSn for the processes (H2O)n þ H2 f (H2O)n(H2) at T = 5 K, T = 100 K, T = 160 K, and T = 220 K. Each process is identified by a point given by the value of n, and the points are connected by straight lines.

clusters and reflect the structure of the underlying potential energy surface. Our calculations have used a variety of Monte Carlo-based methods including parallel tempering, the Kirkwood method, and temperature integration. For the temperature integration we have introduced a quadrature approach that significantly minimizes the number of required temperature evaluations for convergence. For pure water clusters the heat capacity anomalies for all sizes provide valuable insight into the temperature-dependent phase changes of these clusters. These anomalies mark the onset of isomerization, and there is a clear trend from stereoisomer isomerizations in the smaller water clusters to structural isomerizations in the larger water clusters. The isomer of demarcation is the hexamer (n = 6), which exhibits both type of changes. For n < 6, the temperature-dependent heat capacity curves reveal only libration-to-free-rotation (LTFR) isomerizations with no “melting” peaks. For n g 6, all the heat capacity curves show evidence of cluster “melting”. The free energy changes as a function of temperature for prototypical pure water clusters have been shown in Figure 11, and then displayed as a function of cluster size, over the full range of n, in Figure 12. The latter plot exhibits a fine structure that clearly depicts a complex hierarchy of stability. An important 4742


The Journal of Physical Chemistry B contribution that helps unravel the differences in the thermodynamic properties of the cluster processes associated with the free energy changes has been the analysis of the dipole orientations in different geometries. The prefixes cis and trans, which are familiar stereoisomerism terms for describing the orientation of functional groups within molecules, have been applied by analogy to the dipole orientations across atomic centers with hydrogen-bonded interactions. Further, the dipole orientations in the geometric subspaces of certain clusters have been found to be highly repulsive, and significantly influence the observed properties. The results of this analysis are summarized in postulates 1 and 2. These postulates may also be useful tools to understand other cluster systems and may help in the understanding of the relative stabilities of real water clusters. For all cluster sizes examined, under the specified conditions, we find the putative global minimum energy structure for (H2O)nH2 to be qualitatively identical to pure water clusters with an externally adsorbed hydrogen molecule. The classical low-temperature heat capacity behaviors of (H2O)nH2 display anomalies associated with the onset of hindered to free translations and rotations as the hydrogen molecule moves along the exterior of the water core. The structures for (H2O)nH2 explain the lack of features in the Gibbs free energy as a function of cluster size for all temperatures explored in this work. The classical simulations presented in this work have been extended to low temperatures where a classical treatment must be suspect. We have justified our work in these low-temperature regions by questioning how the form of the underlying potential energy surface manifests the computed classical properties. The extent to which the observed features are reflected in real water clusters awaits a fully quantum treatment.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Present Addresses †

Department of Chemistry, The College of The Bahamas, P.O. Box N4912, Nassau, N.P. Bahamas.

’ REFERENCES (1) Dyke, T. R.; Mack, K. M.; Muenter, J. S. J. Chem. Phys. 1977, 66, 498. (2) Xantheas, S. S.; Dunning, T. H. J. Chem. Phys. 1993, 99, 8774. (3) Xantheas, S. S. J. Chem. Phys. 1994, 100, 7523. (4) Pedulla, J.; Kim, K.; Jordan, K. Chem. Phys. Lett. 1998, 291, 78. (5) Pedulla, J.; Jordan, K. Chem. Phys. 1998, 239, 593–601. (6) Rodriguez, J.; Laria, D.; Marceca, E. J.; Estrin, D. A. J. Chem. Phys. 1999, 110, 9039. (7) Nigra, P.; Carignano, M. A.; Kais, S. J. Chem. Phys. 2001, 115, 2621. (8) Tharrington, A. N.; Jordan, K. D. J. Phys. Chem. A 2003, 107, 7380. (9) Kabrede, H.; Hentschke, R. J. Phys. Chem. B 2003, 107, 3914–3920. (10) Shin, S.; Son, W.-J.; Jang, S. J. Mol. Struct.: THEOCHEM 2004, 673, 109–113. (11) Cui, J.; Liu, H.; Jordan, K. D. J. Phys. Chem. B 2006, 110, 18872–18878. (12) McCarthy, V. N.; Jordan, K. D. Chem. Phys. Lett. 2006, 429, 166–168. (13) Langley, S. F.; Curotto, E.; Freeman, D. L.; Doll, J. D. J. Chem. Phys. 2007, 126, 084506–15.

ARTICLE

(14) Frantsuzov, P. A.; Mandelshtam, V. A. J. Chem. Phys. 2008, 128, 094304. (15) Hock, C.; Schmidt, M.; Kuhnen, R.; Bartels, C.; Ma, L.; Haberland, H.; v.Issendorff, B. Phys. Rev. Lett. 2009, 103, 073401. (16) Gelman-Constantin, J.; Carignano, M. A.; Szleifer, I.; Marceca, E. J.; Corti, H. R. J. Chem. Phys. 2010, 133, 024506. (17) Abraham, F. F. Homogeneous Nucleation Theory; Academic Press: New York, 1974. (18) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (19) Hartke, B. Phys. Chem. Chem. Phys. 2003, 5, 275–284. (20) Vitek, A.; Kalus, R.; Paidarova, I. Phys. Chem. Chem. Phys. 2010, 12, 13657. (21) Lynden-Bell, R. M.; Wales, D. J. J. Chem. Phys. 1994, 101, 1460. (22) Predescu, C.; Frantsuzov, P. A.; Mandelshtam, V. A. J. Chem. Phys. 2005, 122, 154305. (23) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300. (24) Egorov, A. V.; Brodskaya, E. N.; Laaksonen, A. J. Chem. Phys. 2003, 118, 6380. (25) Lukyanov, S.; Zidi, Z.; Shevkunov, S. Chem. Phys. 2007, 332, 188–202. (26) Sunden, A. E. K.; Stochkel, K.; Panja, S.; Kadhane, U.; Hvelplund, P.; Nielsen, S. B.; Zettergren, H.; Dynefors, B.; Hansen, K. J. Chem. Phys. 2009, 130, 224308. (27) Douady, J.; Calvo, F.; Spiegelman, F. Euro. Phys. J. D 2009, 52, 47–50. (28) Hartke, B. J. Chem. Phys. 2009, 130, 024905. (29) Alavi, S.; Ripmeester, J. A.; Klug, D. D. J. Chem. Phys. 2005, 123, 024507. (30) Lee, J. K.; Barker, J. A.; Abraham, F. F. J. Chem. Phys. 1973, 58, 3166. (31) Sabo, D.; Freeman, D. L.; Doll, J. J. Chem. Phys. 2005, 122, 094716. (32) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (33) Neirotti, J. P.; Calvo, F.; Freeman, D. L.; Doll, J. D. J. Chem. Phys. 2000, 112, 10340. (34) Neirotti, J. P.; Freeman, D. L.; Doll, J. D. Phys. Rev. E. 2000, 62, 7445. (35) Marinari, E.; Parisi, G. Europhys. Lett. 1992, 19, 451. (36) Geyer, C. J.; Thompson, E. A. J. Am. Stat. Assoc. 1995, 90, 909. (37) Falcioni, M.; Deem, M. J. Chem. Phys. 1999, 110, 1754. (38) CRC Handbook of Chemistry and Physics, 83rd ed.; CRC Press: Boca Raton, FL, 2002. (39) Allen, M.; Tildesley, D. Computer simulations of liquids; Oxford University Press: Oxford, UK, 1987. (40) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269–6271. (41) Zwanzig, R. W. J. Chem. Phys. 1954, 22, 1420. (42) Jarzynski, C. Phys. Rev. Lett. 1997, 78, 2690. (43) Jarzynski, C. Phys. Rev. E. 1997, 56, 5018. (44) Frenkel, D.; Smit, B. Understanding Molecular Simulations; Academic: San Diego, CA, 1996. (45) Boresch, S.; Karplus, M. J. Phys. Chem. A 1999, 103, 103–118. (46) Boresch, S.; Karplus, M. J. Phys. Chem. A 1999, 103, 119–136. (47) Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Chem. Phys. Lett. 1994, 222, 529–539. (48) Reiss, H. Adv. Chem. Phys. 1994, 9, 1–84. (49) Simonson, T. Mol. Phys. 1993, 80, 441–447. (50) Resat, H.; Mezei, M. J. Chem. Phys. 1993, 99, 6052. (51) Freeman, D. L.; Doll, J. D. J. Chem. Phys. 1985, 82, 462. (52) Gautschi, W. Math. Comput. 1968, 22, 251–270. (53) Gautschi, W. Commun. ACM 1968, 11, 432–436. (54) Stillinger, F. H.; Weber, T. A. Phys. Rev. A. 1983, 28, 2408. (55) Kuipers, J. Quaternions and rotation sequences; Princeton University Press: Princeton, NJ, 1999. (56) Rapaport, D. The art of molecular dynamics simulation; Cambridge University Press: Cambridge, UK, 2003. 4743



ARTICLE

(57) Karney, C. J. Mol. Graphics Model. 2007, 25, 595–604. (58) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes, 2nd ed.; Cambridge University Press: Cambridge, UK, 1992. (59) Koonin, S.; Meredith, D. Computational physics; Westview Press: Boulder, CO, 1990. (60) Wales, D. J.; Doye, J. P. K.; Dullweber, A.; Hodges, M.; Naumkin, F.; Calvo, F.; J. Hernandez-Rojas,; Middleton, T. The Cambridge Cluster Database http://wwwwales.ch.cam.ac.uk/CCD.html. (61) Taketsugu, T.; Wales, D. J. Mol. Phys. 2002, 100, 2793–2806. (62) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 1479. (63) Xantheas, S. S.; Burnham, C. J.; Harrison, R. J. J. Chem. Phys. 2002, 116, 1493. (64) Matro, A.; Freeman, D. L.; Topper, R. Q. J. Chem. Phys. 1996, 104, 8690. (65) Fowler, J. E.; Schaefer, H. F. I. J. Am. Chem. Soc. 1995, 117, 446–452. (66) Pinto, Y.; Hel-Or, H.; Avnir, D. J. Chem. Soc., Faraday Trans. 1996, 92, 2523–2527. (67) Xantheas, S. S.; Dunning, T. H. J. Chem. Phys. 1993, 98, 8037. (68) Wales, D. J. Energy Landscapes; Cambridge University Press: Cambridge, UK, 2003. (69) Smith, B. J.; Swanton, D. J.; Pople, J. A.; Schaefer, H. F.; Radom, L. J. Chem. Phys. 1990, 92, 1240. (70) Takahashi, M.; Watanabe, Y.; Taketsugu, T.; Wales, D. J. J. Chem. Phys. 2005, 123, 044302. (71) Holden, G. L. “Free energy changes in water clusters”. Ph.D. Thesis, University of Rhode Island, Kingston, RI, 2010. (72) Tsai, C. J.; Jordan, K. D. J. Chem. Phys. 1993, 99, 6957. (73) Wales, D. J.; Ohmine, I. J. Chem. Phys. 1993, 98, 7245. (74) Asare, E.; Musah, A.; Curotto, E.; Freeman, D. L.; Doll, J. D. J. Chem. Phys. 2009, 131, 184508–11. (75) Hernandez-Rojas, J.; Gonzalez, B. S.; James, T.; Wales, D. J. J. Chem. Phys. 2006, 125, 224302. (76) Vega, C.; Abascal, J. L. F. J. Chem. Phys. 2005, 123, 144504. (77) Freeman, D. L.; Doll, J. D. Adv. Chem. Phys. 1988, 70 (part 2), 139. (78) Stogryn, D. E.; Stogryn, A. P. Mol. Phys. 1966, 11, 371.

4744


Monte Carlo Investigation of the Thermodynamic Properties of (H2O)n

Recommend Documents