The Critical Role of Anharmonicity in Aqueous Ionic Clusters Relevant

J. Phys. Chem. C 2007, 111, 4977-4983

4977

The Critical Role of Anharmonicity in Aqueous Ionic Clusters Relevant to Nucleation Shawn Kathmann,* Gregory Schenter, and Bruce Garrett Chemical Sciences DiVision, Pacific Northwest National Laboratory, Richland, Washington 99352 ReceiVed: NoVember 10, 2006; In Final Form: February 7, 2007

Calculating accurate free energies of aqueous ionic clusters is important in many areas of chemical physics. A common practice used in ab initio and other methods to determine the cluster free energy is to use a single energetically favorable configuration and then assume the harmonic approximation to the potential energy surface. In the present study, we compare and contrast harmonic free energies to those determined from anharmonic configurational sampling and how both compare to experiment. It is found that anharmonicity plays a critical role in determining accurate free energies for the clusters underlying ion-induced nucleation. Consequently, energies, geometries, and frequencies obtained by using a single or multiple configurations within the harmonic approximation lead to inaccurate aqueous cluster thermodynamics.

Introduction Ion-induced nucleation has a long history with several outstanding problems in both theory and experiment.1-31 One important problem is the relative enhancement of the nucleation of water onto positive or negative ions. Laboratory studies of ion-induced nucleation have produced a variety of results and interpretations (sometimes contradictory) depending on the physical conditions and chemical reactions involved. Our previous work concluded that a chemical distinction can explain the origin of discrepancies between experiments.31 The need for theoretical studies has been underscored by Kane et al.,32 who stated that “A complete understanding of the role of various molecular properties in the ion nucleation process is necessary for the development of proper molecular theories of nucleation, with the implicit assumption that these theories will provide a better picture of how ion-molecule and intermolecular interactions lead to the rich dynamic behavior of condensed phase systems.” The focus of this work is to show that proper accounting of anharmonicities in the intermolecular interactions and statistical mechanics plays an important role in quantifying cluster free energies. Furthermore, we will illustrate that the commonly used harmonic approximation is not adequate for the clusters that control nucleation. The goal of molecular theories of nucleation is to obtain nucleation rates starting from the molecular interaction potential.33,34 Dynamical Nucleation Theory (DNT) is a molecularlevel treatment of nucleation kinetics based upon variational transition state theory to obtain rate and equilibrium constants for the underlying clustering reactions. From these quantities, nucleation rates can be obtained as a function of temperature and supersaturation.31,35-41 DNT represents a theoretical framework that connects the interaction potential to the nucleation rate. A molecular theory of nucleation should address three components: (1) accurate descriptions of the interaction potentials, (2) adequate simulation procedures to evaluate cluster thermodynamics and kinetics, and (3) a consistent theoretical approach that connects the interaction potentials to the nucleation rate. Our objective here centers on issues 1 and 2 as they relate to the requirements of computational chemistry and molecular simulation methods for use within DNT. * Author to whom correspondence should be addressed. Tel.: 509-3752870. Fax: 509-375-4381. E-mail: [email protected].

Interaction potentials based on either empirical or ab initio data do a reasonable job of reproducing some properties, but fail for other properties for a number of reasons to be discussed below. However, our recent studies38,40 on the nature of nucleation sensitivity demonstrated that subtle features of the interaction potentials can influence the nucleation rate profoundly and thus interaction potentials should be carefully chosen.42 Empirically derived pair potentials, parametrized to reproduce a particular set of observables, implicitly include many-body effects as well as zero-point energy (ZPE). Quantum derived pair potentials (i.e., those derived from electronic structure calculations), parametrized to ground state binding energies, implicitly include many-body effects, but do not include ZPE. Moreover, quantum derived potentials must somehow account for quantum effects on nuclear motion, such as ZPE effects, either by direct simulation (e.g., by Feynman path integrals) adding to the already large computational cost or by implicit parametrization.42 In the present study we focus on the importance of anharmonicity and avoid the complication of quantum effects on nuclear motion by employing empirically derived interaction potentials for which classical simulations are appropriate. Methods: Theory The presence of ions can lead to the formation of stable aqueous ionic clusters at ambient conditions.43 For small precritical aqueous ionic clusters, the reversible work of formation decreases to a local minima (via successive hydration steps) before rising up to the ion-induced nucleation barrier: the top of the barrier being the location of the critical cluster. For clusters larger than the critical size, the work of formation continues to decrease as the clusters grow into macroscopic droplets. The number of hydration steps required to reach the hydration minima is dictated by the interactions between water and the appropriate ion. In general, this occurs when adjacentsized hydration free energy differences approach those of pure water clusters. The ion-induced nucleation barrier is the difference between the free energy at the hydration minima and the free energy of the critical cluster. For example, clustering reactions around a Na+ ion might yield a minimum for the eighth hydrate whereas the critical cluster may contain 40 water

10.1021/jp067468u CCC: $37.00 © 2007 American Chemical Society Published on Web 03/14/2007

4978 J. Phys. Chem. C, Vol. 111, No. 13, 2007

Kathmann et al.

molecules. Thus, understanding ion-induced nucleation at atmospheric conditions requires being able to calculate the free energies of aqueous ionic clusters containing tens to hundreds of water molecules, not simply the initial ion hydration steps.27,28 The pre-critical clusters governing the ion-induced nucleation barrier (in the example above, those clusters containing between 8 and 40 water molecules) are thermodynamically unstable: they evaporate faster than they grow by condensation. In this study we employ empirically derived interaction potentials that implicitly account for quantum effects on nuclear motion as discussed above. Therefore we use a classical statistical mechanical formulation that allows a direct assessment of the importance of anharmonicity. All the aqueous ionic clusters sample configurations dictated by the ambient temperature and available configuration space defined in the statistical mechanical formalism. For clusters important in nucleation, the statistical mechanics must account for anharmonicity by sampling many cluster configurations. It is important to understand that we are referring to two different types of anharmonicity: (1) local anharmonicity of the vibrations for a given cluster configuration and (2) global anharmonicity resulting from sampling between the large number of configurations available within the relevant volume of configuration space. The i-cluster anharmonic canonical partition function (with center-of-mass translation removed) is given by

qi(V,T) )

γi i!

∫ ‚ ‚ ‚ ∫V dri exp[-U(ri)/kBT]Θ[-S(ri)]

(1)

where γ is the inverse de Broglie wavelength cubed, kB is Boltzmann’s constant, T is the temperature, and U(ri) is the interaction potential (or PES) for a configuration ri. Equation 1 is written for atoms and generalization to internally polyatomic molecules is straightforward. Note that in our previous implementation of DNT the above multidimensional integral in eq 1 is performed over a spherical configuration space denoted by the volume V where the cluster center-of-mass is coincident with the center of the sphere placed at the origin; however, this is not a limitation of DNT as other definitions of the dividing surface can be chosen. The constraining sphere characteristic function Θ is given by i

Θ[-S(ri)] )

θ(rcut - |rj - Ri|) ∏ j)1

(2)

where θ(x) is the Heaviside step function, and Ri is the cluster center-of-mass. The i-cluster anharmonic Helmholtz free energy is

Ai(V,T) ) -kBT ln[qi(V,T)]

(3)

We now turn to how the anharmonic cluster free energies Ai are calculated. The Finite Time Variational External Work44-46 method has been developed into a practical and efficient Monte Carlo procedure38 for calculating cluster Helmholtz free energies. This method is particularly useful since it provides not only cluster Helmholtz free energies but also error estimates that can be systematically improved. In the External Work method, the free energy is determined by an incremental summation of (1) the energy exchanged between the cluster and the heat bath, i.e., changing the populations of different energy states on a fixed PES, and (2) the variation in the cluster Hamiltonian, i.e., leaving the populations fixed but changing the underlying PES. This method is an amalgam of concepts from Clausius, Helmholtz, Metropolis, and Schro¨dinger. The key concept can

Figure 1. The Finite Time Variational External Work method is used to calculate the anharmonic Helmholtz free energy difference (long dashed line) at temperature T between an ideal gas reference state (dotted line U ) 0) and the fully interacting aqueous ionic cluster (solid line U ) U(ri)) at a particular radius, rcut, of the confining sphere.

be traced back to Schro¨dinger, who interpreted doing work on a system as being equivalent to changing the system Hamiltonian. From the First and Second Laws of Thermodynamics it can be shown that a reVersible free energy change is equivalent to the work done on the system. For any other path connecting the thermodynamic states of interest the free energy difference will always be less than the work. By traversing forward and reverse paths between the thermodynamic states of interest, we can obtain an estimate of the error in the free energy difference calculation. Here, the anharmonic free energy difference (see Figure 1) represents the free energy between a convenient reference state (an ideal gas, U(ri) ) 0) and a desired state (a fully interacting system, U ) U(ri)) calculated at some convenient value of the constraining radius rcut. Thus, during the course of the calculation a tally is kept of the changes in the PES (which is gradually turned from ideal, U(ri) ) 0, and interacting, U ) U(ri), states in the forward and reverse directions) and the changes in the heat flux in and out of the cluster. In general, the more Monte Carlo steps that are taken to sample between the two states the closer one is to a fully reversible path. Our simulations were run long enough so that convergence to a reversible path was reached. From these quantities, the cluster free energy can be determined (see Figure 1). For aqueous ionic clusters, we consider i water molecules clustering around a single ion (e.g., Na+ and Cl-). The canonical partition function, QIG, for an ideal gas (U ) 0) cluster of i rigid water (w) molecules plus an ion (X), where the cluster center-of-mass is placed at the origin of a spherical constraining volume of radius rcut, is expressed as IG Q(X‚(H 2O)i 3/2

i γ wV

) C(22i+1πi-1

∫0∞ dy y-(3i+1)

[sin ay - ay cos ay][sin by - by cos by]i) (4) where

Calculating Free Energies of Aqueous Ionic Clusters

C)

[ ]

i γXγi-1 w Ωw MTOT

i3/2i!

J. Phys. Chem. C, Vol. 111, No. 13, 2007 4979

3

mXmiw

where γ is the inverse thermal de Broglie wavelength cubed, V is the total translational volume available to the center-of-mass of the cluster, Ω is the rotational partition function, MTOT is the total cluster mass, a ) mXrcut, b ) mwrcut, m is the mass of the ion or water molecule, and h is Planck’s constant. The aqueous ionic cluster Helmholtz free energy is given as

AXIG(‚(H2O)i(rcut,T)

[ ]

) -kBT ln

QXIG(‚(H2O)i i3/2γwV

(5)

The Helmholtz free energy of the same aqueous ionic cluster with interactions (i.e., with U ) U(ri)), AXINT (‚(H2O) , is calculated i via Monte Carlo simulations, using the Finite Time External Work Method discussed above. The above integral in eq 4 is evaluated numerically for a given number of water molecules, ion and water masses, and constraining radius rcut [the values of rcut used for the aqueous ionic clusters i )1-7 are the following: 2.8, 4.3, 5.8, 5.8, 6.3, 6.8, 7.5 Å]. Convergence of the Helmholtz free energy, AXINT (‚(H2O) , is obtained by a series of i increasingly longer independent Monte Carlo simulations such that the free energy changes by less than 0.1 kcal/mol. The total cluster Helmholtz free energy is obtained by summing the free energies of the ideal gas cluster and the interacting cluster IG INT AXTOT (‚(H2O) ) AX(‚(H2O) + AX(‚(H2O) (see Figure 1). i i i

o Figure 2. Obtaining Gibbs free energies differences ∆Gi,i-1 from the EQ Van’t Hoff plot requires fitting equilibrium constants Ki,i-1 ) Ki,i-1 /Ko measured over a narrow temperature range to a straight line. The slope and intercept are then used to extrapolate the straight line to the temperature of interest, which may lie far beyond were the data were measured.

uncertainty in the resulting free energy difference could easily be a few kilocalories/mole. To transform our results to the standard reference state used in the experiments we start with the expression

Ii Ii-1

) exp[-β(∆Ai,i-1 - µ1)]

where the ambient monomer chemical potential is

[ ]

Results and Discussion Na+(H2O)i and Cl-(H2O)i Clusters. Here we provide a comparison of our calculations on Na+(H2O)i and Cl-(H2O)i to experiments at 243 and 298 K. The purpose of this comparison is threefold: (1) to compare our anharmonic free energy differences to experiment,47,48 (2) to show how the anharmonicity changes as the temperature is increased, and (3) to compare against the harmonic results, using a single or multiple minima. The experimental Gibbs free energies are determined from the Van’t Hoff equation

-

[ ] [

]

o EQ Ki,i-1 Ii ∆Gi,i-1 ) ln ) ) ln o kBT K Ii-1PLKo

o ∆Gi,i-1 ,

EQ Ki,i-1 ,

o ∆Hi,i-1 ,

o o ∆Hi,i-1 1 ∆Si,i-1 + (6) kB T kB o ∆Si,i-1

where and are the Gibbs free energy difference, equilibrium constant, enthalpy difference, and entropy difference between i and i - 1 clusters, respectively, and Ko is the value of the reaction quotient at the standard state. Ii is the experimental ion intensity for an i-cluster. Note that o EQ , Ki,i-1 , the temperature dependence of the quantities ∆Gi,i-1 o o ∆Hi,i-1, ∆Si,i-1, and Ii is implicit. The experimental ambient vapor pressure of the water monomers in the chamber is PL ) kBTN1/V, where N1/V is the ambient water monomer concentration. We note that comparing the computed free energy differences (using electronic structure, empirical potentials, etc.) to those derived from experimental data should be approached with some caution as the experiments are performed over a narrow range of temperatures, typically far away from the temperature of interest (see Figure 2). Thus, extrapolation is used to obtain enthalpies and entropies at the desired temperature. Depending upon how far one has to extrapolate, the

(7)

µ1 ) kBT ln

N1 γΩV

(8)

and β ) 1/kBT. The ∆Ai,i-1 in eq 7 can be obtained from X( ) the calculated Helmholtz free energy differences, ∆Ai,i-1 TOT TOT AX(‚(H2O)i - AX(‚(H2O)i-1, by the following relation that properly accounts for cluster translation:

i 3 X( ∆Ai,i-1 ) ∆Ai,i-1 - kBT ln 2 i-1

[ ]

(9)

Combining eqs 6-9 gives

[

]

o exp(-β∆Ai,i-1 + βµ1) ∆Gi,i-1 ) ln kBT PLKo

Upon simplification,

[ ]

o ) ∆Ai,i-1 - kBT ln ∆Gi,i-1

N1/γΩV P LK o

(10)

)

∆Ai,i-1 + kBT ln[KoγΩkBT] (11)

For example, water monomers at T ) 243 K: γ ) 54.3 Å-3, Ω ) 31.7, kBT ) 0.483 kcal/mol, Ko ) 1 atm-1, which gives o o ∆Gi,i-1 (T)243K) ) ∆Ai,i-1 (T)243K) + 8.63 kcal/mol (12)

By using eq 11 a valid comparison with experiment can now be made. The Lennard-Jones plus Coulomb interaction parameters used in the present study are the same as those employed in our previous study31 on ion-induced nucleation: H2O interaction is described by the TIP4P potential, Na+ ) small positive ion (q ) +1e, σ ) 2.35 Å,  ) 0.13 kcal/mol) and

4980 J. Phys. Chem. C, Vol. 111, No. 13, 2007

Kathmann et al.

Figure 3. Free energy differences between clusters X(H2O)i and X(H2O)i-1 at 298 and 243 K for Na+ (circles and up triangles) and Cl- (squares and down triangles). Anharmonic results (filled symbols with solid lines) and harmonic results (open symbols and dashed lines) are compared with experimental results for Na+ from Dzidic et al. (up triangles with dotted lines) and for Cl- from Hiraoka et al. (down triangles with dotted lines).

Cl- ) large negative ion (q ) -1e, σ ) 4.45 Å,  ) 0.10 kcal/mol).49 Monte Carlo simulations were performed by using the External Work method to calculate the anharmonic Helmholtz free energies. These simulations used 107 configurations per water molecule in the aqueous ionic cluster. Figures 3 and 4 show how our anharmonic free energy differences compare with experiments at 243 and 298 K, respectively. The uncertainty in our anharmonic free energies is (0.2 kcal/mol. We now test the validity of free energy differences calculated using the rigid-rotor harmonic oscillator approximation by comparing them to the full anharmonic free energy differences. Rigid-Rotor Harmonic Oscillator Approximation. The canonical partition function for an i-cluster can be expressed50 as

q ) qtransqrot-vib

[

]

2πMkBT 2

h

3/2

V

1 V

∫ ‚ ‚ ‚ ∫V dri exp[-βU(ri)]χ(V)

qrot )

where U(ri) is the potential energy function governing the motion of the nuclei within the cluster with respect to the cluster center-of-mass, and χ(V) is a function of the confining volume V that restricts the multidimensional integral in eq 15 to relevant regions of the i-cluster configuration space. The relevant cluster configuration space must be defined in such a way as to avoid configurations of nuclei that cannot be characterized as clusters,

[

]

T3 π1/2 σ Θ1Θ2Θ3

1/2

(17)

where σ is the rotational symmetry number of the single cluster configuration, and the Θk (k ) 1, 2, 3) are the principal rotational temperatures given by

Θk )

[ ] h2 8π2IkkB

1/2

(18)

where Ik are the principal moments of inertia of a single cluster configuration around the kth principal axis. The classical vibrational partition function is 6i-6

Qclass vib )

(15)

(16)

The rotational partition function is

(14)

where M is the total mass of the cluster. The coupled rotationalvibrational partition function for a cluster with fixed center-ofmass can be expressed as

qrot-vib )

uncoupled qrot-vib ≈ qrotqvib ) qRRHOA

(13)

where the i-cluster Helmholtz free energy is A(i,V,T) ) - kBT ln[q], and where V is the total volume available to the i-cluster center-of-mass. The translational partition function is obtained by integrating over the i-cluster center-of-mass variables to yield

qtrans )

in other words, a dividing surface must be defined. A general PES of 3i - 3 dimensions (after removing 3 degrees of freedom for overall cluster translation) can be obtained in a variety of ways (e.g., analytic PES and electronic structure calculations). The most accurate energies are obtained from high-level electronic structure methods (e.g., MP2, CCSD(T), and FullCI) with very large basis sets. Density Functional Theory electronic structure methods are becoming commonplace in an attempt to simplify and expedite the energy evaluations. However, given the large number of degrees of freedom encountered in the clusters relevant to nucleation it is computationally intractable to obtain energies for the millions of cluster configurations needed to populate the canonical ensemble from electronic structure calculations. For example, if 10 electronic structure energies are needed for each degree of freedom underlying the PES for a 10 water molecule cluster, this would require 103i-3 ) 1087 calculationssthis is why empirical potentials are so useful in applying statistical mechanics to the clusters involved in nucleation. In light of these computational difficulties, it is common for practitioners of electronic structure to use the rigid-rotor harmonic oscillator approximation (RRHOA) to simplify the multidimensional configurational partition function so that the coupled rotational-vibrational motion of the nuclei becomes separable yielding convenient analytic results valid only for a single cluster configuration. This means that out of the large number of possible liquid-like and dissociative cluster configurations, a single energetically faVorable configuration is chosensa stationary point where all the energy gradients are zero. The resulting uncoupled rotational vibrational partition function is

[kBT/hνj] ∏ j)1

(19)

where νj are the normal mode vibrational frequencies (assuming rigid water molecules) for a single cluster configuration obtained from a Hessian calculation. Note that when electronic structure calculations are used to approximate the potential energy surface the vibrational partition function is computed quantum mechanically. As mentioned above, we are using empirically derived effective potentials and a meaningful comparison of harmonic and anharmonic results is obtained by using classical statistical mechanics for both. A common misperception is that the reasonable accuracy of the RRHOA for predicting spectroscopic properties of molecules validates the RRHOA in general. However, quantities such as

Calculating Free Energies of Aqueous Ionic Clusters

J. Phys. Chem. C, Vol. 111, No. 13, 2007 4981

Figure 4. (a) At low temperatures a single energetic minima may dominate the free energy; however, at higher temperatures the system samples many minima. (b) When more water molecules are included, there are many minima to choose from, and the lowest energy structure may not be the structure with the lowest free energy. For example, at 243 K, the Na+(H2O)3 configuration (c) has the lowest energy (-61.6 kcal/mol) and the lowest Helmholtz free energy (-67.5 kcal/mol). Conversely, the Cl-(H2O)7 configuration (d) does not have the lowest energy (-95.8 kcal/mol) but does have the lowest Helmholtz free energy (-92.8 kcal/mol) compared to the configuration (e) that has the lowest energy (-99.3 kcal/mol) but not the lowest Helmholtz free energy (-90.7 kcal/mol).

vibrational energy splittings between the ground and first excited vibrational states probe only the low-energy portions of the potential. Molecular clusters important in nucleation have a multitude of low-frequency modes, and the RRHOA for spectroscopy of these types of systems has not been established. In addition, the energy range probed by spectroscopic constants such as anharmonic frequencies is much lower that the energies probed in the configurational sampling done to compute the cluster thermodynamic properties required for nucleation. Some improvements to the RRHOA can be accomplished by using perturbation theory to obtain corrections or constants to account for local anharmonicity, centrifugal distortion, and rotationalvibrational coupling. Anharmonic corrections to the RRHOA extend the range of validity of PES but not enough for nucleation. At the conditions relevant to nucleation, nuclei deviate from stationary points on the PESsi.e., the cluster samples liquid-like or quasicrystallite statessrendering the RRHOA and perturbation corrections invalid. As stated previously, thermal averaging over many configurations is necessary and can be calculated directly by evaluating the multidimensional configurational partition function in eq 15. Note that evaluation of eq 15 is very different from eqs 16-19 even if perturbation theory is used to give and first and second order corrections to the RRHOA. Using the same interaction potentials, we can test how well the RRHOA performs for the cluster free energies calculated anharmonically. These results are shown in Figure 3. We compare our classical anharmonic free energy differences to the classical RRHOA free energy differences using frequencies obtained from the global minima. The RRHOA seems to work well for the first few water molecules around the smaller sized Na+ where the binding energy is largest. But the larger sized Cl- ion does not allow the water molecules to bind harmonically. The disagreement between the anharmonic and RRHOA free energies will increase as the number of water molecules increases and/or the temperature is increased, simply because the volume of relevant configuration space becomes larger51,52

(see Figure 3). Furthermore, this situation cannot be remedied by adding more stable conformations to the partition function51 within the RRHOA (sometimes called the “superposition” approximation n

] ∑ QHarm j

Q ) lim [

nf∞ j)1

since this approach neglects important regions of the anharmonic configuration space and will be discussed below. Figure 4a,b shows how the free energy may be dominated by energetic minima for the smallest clusters at low temperatures; however, at higher temperatures or with more water molecules the system samples many higher lying regions of the potential energy surface. As an example of how the energetic minimum for Cl-(H2O)7 may not necessarily be the configuration with the lowest free energy, see Figure 4c,d,e. Our previous work on ion-induced nucleation31 reported anharmonic cluster free energy differences for clusters containing up to 40 water molecules. However, comparison to the RRHOA free energies for clusters containing more water molecules becomes difficult because finding the lowest minima becomes increasingly cumbersomes this is why we stopped the RRHOA analysis at seven water molecules. But, it should be noted that the cluster free energies obtained by using multiple configurations are more negative than free energies obtained by using only a single configurations this feature is not apparent from the differences in cluster free energies. Under what circumstances the RRHOA remains valid is a matter of future research; however, one can say, in general, that it is temperature, system, and size dependent. In the present study, we find that anharmonicity becomes important after just a few water molecules for Na+(H2O)i and immediately for Cl-(H2O)i. We now test how well the superposition RRHOA approximation free energy differences compare with the fully anharmonic results. An ensemble of cluster minima for finding RRHOA configurations was generated by using the following thermal-

4982 J. Phys. Chem. C, Vol. 111, No. 13, 2007

Kathmann et al. Conclusions Prediction of accurate nucleation rates requires (1) valid cluster energetics, (2) valid statistical mechanical methods for obtaining the appropriate cluster thermodynamics (e.g., free energies and equilibrium constants), and (3) a valid chemical kinetics formalism to calculate evaporation and condensation rate constants used in the kinetics mechanism. All three are essential to obtain accurate predictions of absolute nucleation rates. It has been shown that deviations in the underlying PES, whether from ab initio or empirical potentials, can introduce significant errors in nucleation rate prediction. Even though electronic structure calculations can provide accurate energetics for individual cluster configurations, the assumption of the RRHOA on a single or even multiple cluster configurations can invalidate any predictability even when the highest level ab initio calculations are performed. We have shown that anharmonicity (both local and global) in the statistical mechanics of aqueous ionic clusters plays a critical role in the prediction of ion-induced nucleation thermodynamics.

Figure 5. The error due to the difference between the calculated and calc expt experimental ∆Gi,i-1 (error ) ∆Gi,i-1 - ∆Gi,i-1 ) for Na+ and Cl- at single (squares with 243 and 298 K: single configuration RRHOA ∆Gi,i-1 multiple dashed lines), multiple configuration superposition RRHOA ∆Gi,i-1 anh (diamonds with dotted lines), and the fully anharmonic ∆Gi,i-1 (circles with solid lines).

ization and annealing schedule: A constant energy simulation corresponding to an initial thermalization of velocities at 243 K was used to generate configurations. At 5 ps intervals a configuration was taken as the initial condition for an annealing procedure in which constant energy dynamics was performed and the velocities were set to zero each time the kinetic energy of the system reached a maximum. The cluster configuration quickly converged to a local energetic minimum where a harmonic normal-mode analysis was performed (the gradients were checked to be sufficiently small). One hundred such configurations were generated for the Na+(H2O)i system with i < 5 and the Cl-(H2O)i system with i < 3. One thousand configurations were generated for the Na+(H2O)i system with i > 4 and the Cl-(H2O)i system with i > 2. We show in Figure 5 the resulting errors due to differences between the calculated (single, multiple, fully anharmonic) and experimental free energy differences at 243 and 298 K. The errors in the single and multiple configuration free energy differences range from less than 1 kcal/mol up to 4 kcal/mol. The fully anharmonic free energy differences have errors of about 1.3 kcal/mol or smaller and thus can be considered within chemical accuracy. We have also run explicit molecular dynamics for these clusters at 243 and 298 K. The trajectories resulting from the molecular dynamics simulations are consistent with the breakdown of the RRHOAsthat cluster configurations fluctuate quite far from their energetic minima as more water molecules are added or when the temperature is raised from 243 to 298 K. Performing optimizations from multiple initial geometries to a single energetically favorable minimum may be an insufficient procedure to sample the relevant configuration space. Thus, ab initio or other methods that use a single or multiple energetic minima to calculate RRHOA free energies to deduce cluster populations or other cluster properties at ambient conditions relevant to the atmosphere are questionable and should be checked against anharmonic sampling of the cluster configuration space.

Acknowledgment. This work was supported by the Chemical Sciences Division, Office of Basic Energy Sciences, Department of Energy. The Pacific Northwest National Laboratory is operated by Battelle for the U.S. Department of Energy. References and Notes (1) Wilson, C. T. R. Trans. R. Soc. (London) 1897, A189, 265. (2) Wilson, C. T. R. Trans. R. Soc. (London) 1900, A193, 289. (3) Laby, T. H. Philos. Trans. R. Soc. London 1908, 208, 445. (4) Loeb, L. B.; Kipp, A. F.; Einarson, A. W. J. Chem. Phys. 1938, 6, 264. (5) Scharrer, L. Ann. Phys. 1938, 35, 619. (6) Donnan, F. G. Philos. Mag. 1900, 3, 305. (7) White, D. R.; Kassner, J. L. J. Aerosol Sci. 1971, 2, 201. (8) Kassamali, A. Master’s Thesis, Clarkson College of Technology, 1974. (9) Rabeony, H.; Mirabel, P. J. Phys. Chem. 1987, 91, 1815. (10) Katz, J. L.; Fisk, J. A.; Chakarov, V. M. J. Chem. Phys. 1994, 101, 2309. (11) Katz, J. L.; Lihavainen, H.; Rudek, M. M.; Salter, B. C. J. Chem. Phys. 2000, 112, 8363. (12) Castleman, A. W. J.; Holland, P. M.; Keese, R. G. J. Chem. Phys. 1978, 68, 1760. (13) Castleman, A. W. J.; Keese, R. G. Science 1988, 241, 36. (14) Castleman, A. W. J.; Tang, I. N. J. Chem. Phys. 1972, 57, 3629. (15) Kane, D.; Daly, G. M.; El-Shall, M. S. J. Phys. Chem. 1995, 99, 7867. (16) Thomson, J. J. Conduction of Electricity Through Gases; Cambridge University Press: Cambridge, UK, 1906; Vol. 1. (17) Thomson, J. J. Conduction of Electricity Through Gases, 3rd ed.; Cambridge University Press: Cambridge, UK, 1936; Vol. 1. (18) Zettlemoyer, A. C. Nucleation; Marcel Dekker: New York, 1969. (19) Abraham, F. F. J. Chem. Phys. 1969, 51, 1632. (20) Abraham, F. F. Homogeneous Nucleation Theory; Academic Press: New York, 1974. (21) Yu, F.; Turco, R. P. Geophys. Res. Lett. 2000, 27, 883. (22) Rusanov, A. I.; Kuni, F. M. J. Colloid Interface Sci. 1984, 100, 264. (23) Kusaka, I.; Wang, Z.; Seinfeld, J. J. Chem. Phys. 1995, 103, 89939009. (24) Kusaka, I.; Wang, Z.; Seinfeld, J. J. Chem. Phys. 1995, 102, 913924. (25) Oh, K. J.; Gao, G. T.; Zeng, X. C. Phys. ReV. Lett. 2001, 86, 5080. (26) Brodskaya, E.; Lyubartsev, A. P.; Laaksonen, A. J. Chem. Phys. 2002, 116, 7879. (27) Nadykto, A. B.; Yu, F. Phys. ReV. Lett. 2004, 93, 016101. (28) Nadykto, A. B.; Al, Natsheh, A.; Mikkelsen, K. V.; Yu, F.; Ruuskanen, J. Phys. ReV. Lett. 2006, 96, 125701. (29) Adachi, M.; Okuyama, K.; Seinfeld, J. J. Aerosol Sci. 1992, 23, 327-337. (30) Okuyama, K.; Adachi, M.; Shinagawa, H.; Shi, G.; Seinfeld, J. H. J. Aerosol Sci. 1991, 22, S85. (31) Kathmann, S. M.; Schenter, G. K.; Garrett, B. C. Phys. ReV. Lett. 2005, 94, 116104.

Calculating Free Energies of Aqueous Ionic Clusters (32) Kane, D.; Rusyniak, M.; Fisenko, S. P.; El-Shall, M. S. J. Phys. Chem. A 2000, 104, 4912. (33) Oxtoby, D. W. J. Phys.: Condens. Matter 1992, 4, 7627-7650. (34) Reiss, H. In 15th International Conference on Nucleation and Atmospheric Aerosols; Hale, B. N., Kulmala, M., Eds.; AIP: Rolla, MO, 2000; Vol. 534, p 181. (35) Schenter, G. K.; Kathmann, S. M.; Garrett, B. C. Phys. ReV. Lett. 1999, 82, 3484-3487. (36) Schenter, G. K.; Kathmann, S. M.; Garrett, B. C. J. Chem. Phys. 1999, 110, 7951-7959. (37) Kathmann, S. M.; Schenter, G. K.; Garrett, B. C. J. Chem. Phys. 1999, 111, 4688-4697. (38) Kathmann, S. M.; Schenter, G. K.; Garrett, B. C. J. Chem. Phys. 2002, 116, 5046. (39) Schenter, G. K.; Kathmann, S. M.; Garrett, B. C. J. Chem. Phys. 2002, 116, 4275. (40) Kathmann, S. M.; Schenter, G. K.; Garrett, B. C. J. Chem. Phys. 2004, 120, 9133. (41) Kathmann, S. M. Theor. Chem. Acc. 2006, 116, 169.

J. Phys. Chem. C, Vol. 111, No. 13, 2007 4983 (42) Schenter, G. K. J. Chem. Phys. 2002, 117, 6573. (43) Seinfeld, J. H.; Pandis, S. N. Atmospheric Chemistry and Physics: From Air Pollution to Climate Change; John Wiley & Sons, Inc.: New York, 1998. (44) Reinhardt, W. P.; Hunter, J. E., III J. Chem. Phys. 1992, 97, 15991601. (45) Hunter, J. E. I.; Reinhardt, W. P.; Davis, T. F. J. Chem. Phys. 1993, 99, 6856-6864. (46) Hunter, J. E. I.; Reinhardt, W. P. J. Chem. Phys. 1995, 103, 86278637. (47) Dzidic, I.; Kebarle, P. J. Phys. Chem. 1970, 74, 1466. (48) Hiraoka, K.; Mizuse, S. Chem. Phys. 1987, 118, 457. (49) Smith, D. E.; Dang, L. X. J. Chem. Phys. 1994, 100, 3757. (50) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (51) Brooks, C. L., III; Onuchic, J. N.; Wales, D. J. Science 2001, 293, 612. (52) Doye, J. P. K.; Miller, M. A.; Wales, D. J. J. Chem. Phys. 1999, 110, 6896.