The Thermodynamics and Structure of Hydrated Halide and Alkali Ions BY FARID F. ABRAHAM, MICHAEL R. MRUZIK IBM Research Laboratory, San Jose, California 95193 AND G. MARSHALL POUND Materials Science Department, Stanford University, Stanford, California 94303, U.S.A. Received 12th December, 1975 Gibbs free energies were calculated for the gas phase reaction: ion(H20)~- 1 +H&apour) = ion(H2Oh for the Li+, Na', K+, C1-, and F- ions and for N = 1 to 6. The Monte Carlo method "as used to evaluate the appropriate classical expressions of statistical mechanics by employing the inter- molecular potential functions recently developed from ab initia Hartree-Fock calculations. Enthal- pies and structural information were also calculated. Agreement with experiment is sufficiently good to demonstrate the feasibility of this approach.1. INTRODUCTION Adapting the intermolecular water-water and water-ion potential functions developed from ab initio Hartree-Fock calculations, we demonstrate a successful application of the Monte Carlo method of classical statistical mechanics to evaluate the thermodynamic and structural properties of hydrated halide and alkali ions. This study deviates from heuristic theories of gas phase hydration of ions where the physical cluster is pictured as a " liquid droplet," a " microcrystallite," or some other model construct. While these models provide an insight into the origin of certain features of a physical cluster, they cannot serve to elucidate the cluster's true molecular structure, since this is assumed in one way or another.With the advent of large scale computers, pursuing a first-principles approach for many-body systems has begun where only a form of the intermolecular potential function between two molecules is a s s ~ m e d . ~ Recent examples of this approach have been the molecular dynamics and Monte Carlo simulations of Lennard-Jones atomic clusters, and Monte Carlo simulation of water clusters.8 In Section 2, the Hartree-Fock potential functions for water-water and water-ion interactions are presented. The Monte Carlo method is outlined in Section 3, and expressions for various thermodynamic state variables are derived in terms of quanti- tives that can be measured using the Monte Carlo method. A brief description of recent experiments 9* lo on gas phase hydration of alkali and halide ions is given in Section 4.In Section 5, we present recent theoretical free energies of formation of ion-water clusters for comparison with experiment." In particular, Gibbs free energies have been calculated for the gas phase reaction: i0n(H,0),,~- I-/- H20~vapour) ion(H,O), (1)F A R I D F . ABRAHAM, MICHAEL R . MRUZIK A N D G . MARSHALL P O U N D 35 for the Li+, Na+, K f , C1-, and F- ions and for N = 1 to 6 . Enthalpies, and structural information are also presented. Agreement with experiment is sufficiently good to demonstrate the feasibility of this theoretical approach. 2. WATER-WATER AND WATER-ION POTENTIAL FUNCTIONS Water-water potential functions have been produced which are based on first principles ca1culations.l These functions obtain their form from various point charge models and are adjusted to agree with the ab ivlitio calculations of the Hartree-Fock potential energy surface. The C-XI1 function is the most recent in a series of suc- +Q *H2 FIG.1.-The Beriial and Fowler charge model l2 is employed for the C-XI1 water-water potential function.' The symbol H represents hydrogen atoms; 0 represents oxygen; and M is a fictitious point along the line of symmetry of the water molecule. The OM distance is 0.225 954 A and the OH distance is 0.957 A. The HOH angle is 105". cessively improved formulas of this type. This function uses the Bernal and Fowler charge model l2 illustrated in fig. 1 and is given by: UC-XII = UHartree-Fock+ Udispersion UHartree-Fock = Q2(1/r13+1/r14+ l/r23$.1/r24)+4Q2/r78 -2Q2(1/r18+ 1/r28+ 1/r37+ 1/r47) +a, exp(--b,r,,) (2) + a2[exp(- b2r13) + exP(--b2r14) +exp(- b2r23) + exP(- b2324)I +a3 [exp(-b3r16) + exp(-b,r26) +exp(--b3r35) + exp(- b3r45)1 Udispersion = cl/r566-c2/r586+c3/r4~ (3) where the constants are given by, Q2 = 139.272 kcal A/mol a, = 71533.4 kcal/mol b, = 3.96994/A c1 = 922.781 kcal A6/mol a2 = 779.885 kcal/mol b2 = 3.12544/A c2 = 17283.5 kcal A8/mol a, = 4084.02 kcal/mol b3 = 3.91443/A c3 = 24119.7 kcal bilo/mol The Hartree-Fock energy surface has been modified to include correlation effects which decrease, by making more negative, the total water-water binding energy by 15% to 30%. The dispersion terms were developed from the data of the Quantum Chemistry Group of the University of Warsaw, directed by Prof.W. K0l0s.l~ A rigid geometry is assumed for the water molecule and multibody effects are ignored. Errors in the water-water potential function are not highly significant in this study, since the dominant force in ion-water clusters is the ion-water interaction. For example, the empirical ST2 water-water potential function,14 which has been " effec- tively " adjusted so as to include the average multibody and other effects present in bulk water, gives essentially the same results in this study as does the C-XI1 function." As in the case of the water dimer, water-ion potential functions have been obtained36 THERMODYNAMICS AND STRUCTURE OF HALIDE AND ALKALI IONS from ab initio Hartree-Fock calculations.2 The functions which approximate this Hartree-Fock potential energy surface were based on two arbitrary point charge models: the " simple " ion-water model, which is a modification of the Bernal and Fowler version, and the slightly more complicated " accurate " model, as illustrated in fig.2. The corresponding functions are : / a ) Simple Model Ion PI '10 Ql ( b l Accurate Model FIG. 2.--Charge models employed for the water-ion potential functions.' Point charges are placed on the hydrogen atoms (HI and H2), oxygen atoms (0), along the water molecule line of symmetry (M), and along the bisectors of the OM distance (MI and M2). In (a) all lines are in the plane of the paper, and in (6) all dotted lines are parallel to their respective three dimensional axes. where the numerical values of the constants are listed in tables V and VI of ref, (2). The correlation terms are negligible.We have used the " simple " ion-water poten- tial function for the Li+ and Na+ ions and the " accurate " potential function for F-, Cl-, and K+. For a given ion-water configuration, both functions give almost identical values of potential energy, but the " accurate " version requires slightly more computing time. Although both ion-water potential energy functions were very close to the Hartree-FARID F . ABRAHAM, MICHAEL R . MRUZIK A N D G . MARSHALL P O U N D 37 Fock limit, errors were introduced by assuming a rigid geometry for the water mole- cule (rigid rotor assumption) and by neglecting correlation and multibody effects. The net result of these assumptions and omissions is that the ion-water potential functions should yield values of potential energies which are too positive by about 5% for clusters containing three or less water molecules and which are too negative by as much as 10% for larger clusters containing up to six water molecu1es.l1 To further investigate the effects of these errors, the enthalpy and free energy of reaction (1) for the Na+ ion and N = 6 were determined for a 10% positive change in ion- water and water-water potential functions.The free energy was reduced by 13% to -5.025 kcal/mol, and the enthalpy was reduced by 16% to -10.71 kcal/mol. These estimates are representative of the expected error in the Monte Carlo calculations. 3. THE MONTE CARL0 METHOD In the Monte Carlo calculations, the ion-water clusters are defined by the criterion that all water molecules of the cluster should lie within a spherical volume, V,, centred on the ion.This constraining volume was chosen large enough so that the strong ion-water binding naturally restricted the water molecules to configurations well away from the constraining boundary. Calculations performed with respect to this frame of reference with the ion fixed at the origin were later corrected for the loss of three degrees of freedom associated with cluster translation. This cluster definition did not otherwise affect the canonical averaging process since those excluded con- figurations represented very high energy states, each with a negligible probability of occurrence. The standard state for all calculations was taken as 298 K and 1 atmos- phere.Since Metropolis et aZ.15 first applied the “ Monte Carlo ’’ method to statistical mechanics, this method has been discussed in detail by various author^.^ We will only present a brief review. The method is based on a stochastic process which generates a Boltzmann-weighted chain of configurations of a given N-particle system. The mean value of any function of the system’s coordinates over all configurations in the chain provides an estimate of the canonical ensemble average of that function. For example, the mean potential energy leads to an estimate of the internal energy. The procedure is as follows (see fig. 3). We consider N water molecules in an initial ice cluster configuration of uniform density equal to the experimental water density at temperature T.The cluster is centred on the ion in a constraining spherical boundary with radius R, and volume V, = SV,, where V,, is the volume for bulk liquid water of N molecules at temperature T. Within the framework of our definition of the ion-water cluster, we generate configurations in the following manner: (i) select a molecule at random, (ii) select displacements Ax, Ay, Az, each uniformly distributed on (-A/2, A/2), (iii) select x , y , or z axis at random, (iv) select an angle 8 uniformly distributed on (-w,w); (v) if these displacements move the centre of mass of any water molecule outside the spherical constraining boundary with origin corresponding to the position of the ion, reject the try and accept the old configuration; otherwise (vi) calculate the change in potential energy SUN on displacing the chosen molecule by (Ax, Ay, Az) and rotating it through 8 about the chosen axis; (vii) if SUN is nega- tive, accept the new configuration; otherwise (viii) select a number h uniformly distri- buted on (O,l), (ix) if exp(--GU,/kT) < h, accept the old Configuration; otherwise, (x) the new configuration and the new potential energy become the “current” properties of the system.These rules ensure that averaging over long chains approaches classical canonical averaging, with weighting of configurations proportional to exp(- U,/kT). The38 THERMODYNAMICS AND STRUCTURE OF HALIDE A N D ALKALI IONS The Monte Carlo Procedure 1. Given (old) Configuration Mo lecu I e Moves Out of v 2. Select 3. Displace and Rotate Water Molecule 4.Calculate Energy Change 5. Metropolis Test M.:h t i - random no I 6. Accepted (New) I FIG. 3.-The Monte Carlo algorithm for calculating canonical averages of state variables for ion- water clusters. angle ly was fixed at 0.2 radians and A was chosen so that approximately half the attempted moves were actually made. The change in Helmholtz free energy AFN-l,N for reaction (1) is given by 1 A F N - 1 . N = {l <u,,o(4> +kT In (N- Iy, PdV, (6) >Vc Y(l atm.) where { UHzo(;l)> is the canonical average of the potential energy of the water molecule (monomer) which is being added to the cluster in the following way : the dimensionless coupling parameter, 0 < 3, < 1, scales the interaction of the monomer with the ion and the other N- 1 water molecules.When the interaction is turned on, A == 1, and the monomer is indistinguishable from the other N - 1 water molecules ; but when the interaction is turned off, 3, = 0, and the monomer is completely non-interacting and behaves as an ideal gas. The work of expanding the constraining volume until the unattached monomer (3, = 0) is at a pressure of one atmosphere is described by the pdV term of eqn (6). To compute the canonical average of the monomer potential energy, ( UH,o(A)), 100000 configurations were generated for a number of discrete values of 3,. A numerical integration of these values as a function of 3, furnished the corresponding free energy term. Standard deviations of the thermodynamic variables were obtained from a series of mean values, each of which represented the average of a block con- taining 10 000 configurations.To compensate for the numerical integration difficulties caused by the rapidF A R I D F . A B R A H A M , MICHAEL K. MRUZIK A N D G . MARSHALL POUND 39 increase of (UHzo) as R goes to zero, as illustrated in fig. 4, the following transforma- tion was employed: The value of Z was fairly insensitive to a choice of rn and varied only about 2% for 0.10 < m < 0.75; m was arbitrarily set equal to 0.25 since the average energy for a particle in contact with a fixed energy source is proportional to %-3/4 when the inter- action is dominated by a repulsive Lennard-Jones term of order 2/r.12 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 FIG. 4.--The canonical ensemble average potential energy of a water molecule (UH~O(~)) as a fiinc- tion of the 1 interaction parameter.For the purpose of comparing the theoretical free energies with experiment, the Helmholtz free energy was converted to Gibbs free energy A G N - 1 . N using the relation: A G N - 1 , N = AFN-i,N+A(PV) = AFN-i,N-kT, (8) for the standard state (298 K and 1 atmosphere). The formula used to calculate the changes in internal energy and enthalpy are : A E N - 1 , N = { uN(L=l))-{ u N - 1 ( ] * = I)), (9) AHN-1,N = AEN--,N-kT. (10) 4. EXPERIMENTAL Recent experimental advances 9 9 lo have provided the first empirical thermodynamic data on small ion-water clusters. To make an effective comparison with theoretical calculations, it is appropriate to review the general nature and accuracy of these experiments. Equilibrium constants for the incremental hydration of the alkali metal and halogen ions as per reaction40 THERMODYNAMICS AND STRUCTURE OF HALIDE AND ALKALI IONS (1) were determined as follows (see fig.5 ) : (i) the appropriate gas-phase ions were produced and electrostatically directed to a reaction chamber which contained a known pressure, pH20, (usually about 1 Torr) of water vapour; (ii) a high pressure mass spectrometer recorded the relative concentrations of ion-water clusters as a function of mass number; and (iii) the equilibrium constant, K, and Gibbs free energy, h G N - i , N , were obtained from, Production of Negative Ions Production of Positive Ions (0, c NF3 7 F- i NF:, + 02) 0, Carrier Gas with NF3 I Platinum Grid +- 0 h, Electron Beam or a - particle Source /' __ Anions ~ Focusing Grid Water Molecule To Mass Spectrometer Ion (H201N-, + H,O = Ion (H,O)N FIG.5.-Schematic of experiments for obtaining equilibrium constants for the incremental hydration of the alkali metal and halogen ion as per reaction (l), [ref. (9) and (lo)]. where IN is the concentration of clusters comprised of an ion and N water molecules. Meas- urements of the equilibrium constant as a function of temperature gave the enthalpy of the reaction (l), AHN-i,N, from the van't Hoff relation, Although errors were not estimated in the experimental report, it has been suggested elsewhere that standard deviations of & 1 to i-2 kcal/mol for the enthalpy would not be unreasonable. These estimates presumably include errors such as the dissociation of clusters in vacuum immediately after entering the mass spectrometer.Since the time re- quired for the initial electrostatic acceleration in the mass spectrometer is the same order of magnitude (microseconds) as the period for evaporation of water molecules from the cluster, this effect would decrease the apparent concentration of larger clusters, resulting in more positive free energy measurements. Considerations such as the attainment of equilibrium among the assorted cluster sizes in the reaction chamber and the prevention of anomalous cluster growth, caused by cooling on adiabatic expansion into the mass spectrometer, were examined and probably contributed little to the total error. In addition, different experimental conditions for positive and negative ions complicated a comparison of the error trends of the two groups.Negative ions were produced by ionizing the oxygen to 0; which collided with gaseous NF, or CC14 to create the desired anion. Positive ions were produced by thermionic emission from a platinum filament coated with the appropriate salt. The resulting temperature gradient, which for Li+ began at 900 "CF A R I D F. ABRAHAM, M I C H A E L R . M R U Z I K A N D G. MARSHALL POUND 41 at the filament, was difficult to isolate and would have the effect of making the observed free energies too positive for the positive ions. Cluster sizes larger than six water molecules ( N = 6) could not be studied experimentally due to condensation on the chamber walls which eventually closed the port to the mass spectrometer. A Kf ton (Experimental) A K+ Ion (Monte Carlo) 9 Na' Ion (Experimental) 0 Na' Ion (Monte Carlo) Lit Ion (Experimental) 0 Lit Ion (Monte Carlo) 1 2 3 4 5 6 N, Number of Water Molecules FIG.6.-Enthalpies are given for the reaction M+(H20)N-1+H20 = M+(H~O)N for the cations indicated. All values are for 1 atmosphere and 298 K. L n 5 -24 Q -32 0 CQ- Ion (Experimental) 0 CQ- Ion (Monte Carlo) F' ion (Experimental) F' Ion (Monte Carlo) 1 1 1 1 I 1 1 2 3 4 5 6 #, Number of Water Molecules Fro. 7.-Enthalpies are given for the reaction X - ( H 2 0 ) ~ - ~ + H Z 0 = X-(H20)N for the anions indi- cated. All values are for 1 atmosphere and 298 K.42 THERMODYNAMICS A N D STRUCTURE OF HALIDE A N D ALKALI IONS 5.RESULTS AND DISCUSSION Monte Carlo and experimental lo enthalpies for ion-water clusters containing the F-, Cl-, K+, Na+, Li+ ions and from 1 to 6 water molecules are presented in fig. 6 and 7. As the clusters increase in size, the enthalpy change AHN--l,N for the addition of a water molecule becomes more positive. This is caused by a decrease in ion-water attraction and, in some cases (notably Li+ and Na+), an increase in water-water repulsion with increasing cluster size as illustrated in table 1. The TABLE MONTE CARLO VALUES OF THE AVERAGE ION-WATER POTENTIAL PER WATER MOLE- WATER MOLECULE, ( Uww)/N. ALL VALUES ARE EXPRESSED IN KCAL/MOL CULE, (ulw)/N, IS COMPARED WITH THE AVERAGE WATER-WATER POTENTIAL ENERGY PER 1 -21.89 0 - 10.93 0 -33.44 0 2 - 21.76 0.65 - 10.73 0.08 -33.33 1.25 3 - 21.44 1.39 - 10.58 0.1 1 -33.15 2.74 4 - 21.04 2.08 - 10.08 -0.19 - 32.71 4.15 5 - 19.46 1.71 -9.44 -0.53 -31.94 5.35 6 - 17.98 1.35 -9.18 - 0.74 - 30.76 5.99 1 - 24.32 0 - 16.92 0 2 - 24.24 0.79 - 16.66 0.37 3 - 24.02 1.82 - 16.58 1.10 4 -23.65 2.74 - 16.02 1.39 5 -22.82 3.24 - 14.81 0.98 6 -21.53 3.18 - 13.75 0.57 effect of the high surface-to-volume ratio causes the enthalpy change to become more positive than that for the condensation of a water molecule from the vapour to the bulk liquid (-10.52 kcal/mol) in the cases of K+ and C1-.The Monte Carlo enthalpies are more negative than those from experiment 9* except for the case of C1- where the opposite is true. For the positive ions, discrepancies between the Monte Carlo and experimental data increase in going from K+ to Li+.One notes that the discrepancy between Monte Carlo and experimental enthalpies is generally greater than the estimated error of the Monte Carlo calculations in the cases of Li+ and F-. A comparison of the Monte Carlo free energies with those of experiment 9* lo for the incremental hydration of the different ions is given in fig. 8, 9, and 10. For the positive ions, the Monte Carlo free energies are more negative than experiment, while for the negative ions, the relationship is mixed. The negative deviations of the Monte Carlo from the experimental free energy data for cation clusters increase in passing from Kf to Li+. The experimental free energy for the F- ion displays an anomalous decrease at N = 5. Besides being physically unreasonable, this type of behaviour is not otherwise observed and is most likely an artifact of the experiment.16 For C1-,F A R I D F .ABRAHAM, MICHAEL R . MRUZIK A N D G. MARSHALL POUND 43 which has the best agreement with experiment, the Monte Carlo free energies are slightly more positive than experiment for small clusters (N < 3), and slightly more negative for larger clusters. For the largest clusters studied, N = 6, the experimental free energies approximately converge to the same value (fortuitously close to the value for condensation of vapour to bulk liquid, -2.055 kcal/mol) while the Monte Carlo E -9 CI -12 L 0 u) 0 r ? z u I Ion A KS ton (Experimental) A K+ Ion (Monte Carlo) I I I I 1 1 2 3 4 5 6 ti X , Number of Water Molecules -3 2 .- n -15 0 b) Sodium ton 3 0 Na+ ton (Experimental) 3 Na’ Ion (Monte Carlo) 0 l - 1 I I I I 1 1 2 3 4 5 6 N , Number of Water Molecules FIG.8.-Gibbs free energies are given for the reaction M+(H20)N-l+H20 = Mf(H20), for (a) M+ = K+ and (b) M+ = Na+. All values are for 1 atmosphere and 298 K. free energies are still dependent on the type of ion. One notes that the discrepancy between Monte Carlo and experimental free energies is generally greater than the estimated error of the Monte Carlo calculations in all cases except CI- . In fig. 11 and 12, radial density distributions for oxygen and hydrogen atoms show the formation of a well defined first hydration shell in which the water molecules are constrained to orientations closely resembling the optimum configurations of fig.13. For the anions, there is a mixing of configurations 13(b) and 13(c). Clusters with ions of the same polarity had nearly identical (average) structures differing principally in the distance between the ion and the first hydration shell as determined44 THERMODYNAMICS AND STRUCTURE OF HALIDE AND ALKALI IONS -3 r -30 FIG. 0- I I I I I I - 9.-Gibbs - 21 2 4 1 -27 E Li' Ion (Experimental) Lit Ion (Monte Carlo) Li + ( HzO)~. m I 0 C!2- Ion (Monte Carlo) W F- Ion (Experimental) El F- Ion (Monte Carlo) q- 1 2 c) 5 6 #, Number of Water Molecules FIG. 10.-Gibbs free energies are given for the reaction X-(HZO)N-l+HZO = X-(HzO)N X- = F- and Cl-. All values are for 1 atmosphere and 298 K. All forF A R I D F. A B R A H A M , MICHAEL R . MRUZIK AND G .MARSHALL POUND 45 18 16 for K " ( H ~ o ) , L 4 2 0 0 1 2 3 4 5 6 b) Radial Densities for Na' (H20l6 I,/ Oxwen R, Radial Distance from Icn ( A ) FIG. 1 1 .-Radial densities of oxygen and hydrogen atoms are given for cation clusters of (a) K+(H20)6, (6) Na+(H20)6, and (c) Li+(H20)6.46 THERMODYNAMICS AND STRUCTURE OF HALIDE AND ALKALI IONS 10 9 - 8 - I v1 *% . 7 e 5 - 3 2 6 - L .- .- Q >. 4 - 3 - 2 - 1 - - Ln .- I 1 - b) RadiaCDensities for ce- (H2016 - I a) Radial Densities for F-(H20)6 .(Hydrogen 2 . 3 4 5 6 'R, Radial Distance froin Ion (A) by the first peak in the radial density distribution. Detailed observations of many of the individual configurations generated by the Monte Carlo method indicated no structures reminiscent of the ices or clathrates.In this paper we have performed the first theoretical cluster free energy calculations where a direct verification with experiment is possible and thereby have demonstrated the feasibility of applying the Monte Carlo method in calculating free energies for clusters containing only a few molecules. a) Optimum Cation- C J HB Configuration Water Configuration &!ion l \ l d k r Aiiiuii W3lt.r Aillor1 Wdier - FIG. 13.-Different configurations for ion-water pairs; in (a) and (6) hydrogen atoms are equidistant from the ion and in (c) one hydrogen atom is on a direct line joining the anion and the oxygen atom. All atoms lie in the plane of the paper.FARID F. ABRAHAM, MICHAEL R . MRUZIK A N D G . MARSHALL POUND 47 G . C. Lie and E. Clementi, J . Chettz. Phys., 1975, 62, 2195. W. W. Wood, Physics of Simple Liquids, ed. J. S. Rowlinson, G. S. Rushbrooke and H. N. V. Temperley (North-Holland, Amsterdam, 1968). M. Volmer, Kinetik der Phasenbildung (Theodor Steinlopff Verlag, Dresden, Germany, 1939). E. F. O’Brien and G. W. Robinson, J . Chem. Phys., 1974,61,1050. D. J. McGinty, J . Chem. Phys., 1973, 58, 4733. F. F. Abraham, J. Chem. Phys., 1974, 61, 1221. I. Dzidic and P. Kebarle, J. Phys. Chem., 1970, 74, 1466. M. Mruzik, F. F. Abraham, D. E. Schreiber and G. M. Pound, J. Chem. Phys., to be published. J. D. Bernal and R. H. Fowler, J. Chem. Plzys., 1933, 1, 515. l3 Unpublished data reported in ref. (1). l4 F. H. Stillinger and A. Rahman, J. Chem. Phjs., 1974, 60, 1545. ’ H. Kistenmacher, H. Popkie and E. Clementi, J. Chem. Phys., 1973, 59, 5542. ’ J. K. Lee, J. A. Barker and F. F. Abraham, J. Chern. Phys., 1973, 58, 3166. lo M. Arshadi, R. Yamdagni and P. Kebarle, J. Phys. Chem., 1970,74, 1475. N. Metroplis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. P/?vs., 1953,21, 1087. l6 P. Kebarle, personal communication.