Structure and Spectroscopy of Hydrated Sodium Ions at Different

Article pubs.acs.org/JCTC

Structure and Spectroscopy of Hydrated Sodium Ions at Different Temperatures and the Cluster Stability Rules Jean Jules Fifen†,‡ and Noam Agmon*,† †

The Fritz Haber Research Center, Institute of Chemistry, The Hebrew University of Jerusalem, Jerusalem 91904, Israel Department of Physics, Faculty of Science, The University of Ngaoundere, P.O. Box 454, Ngaoundere, Cameroon

‡

S Supporting Information *

ABSTRACT: The sodium cation plays an important role in several physiological processes. Understanding its solvation may help understanding ion selectivity in sodium channels that are pivotal for nerve impulses. This paper presents a thorough investigation of over 75 isomers of gas-phase Na+(H2O)n=1−8 clusters, whose optimized structures, energies, and (harmonic) vibrational frequencies were computed quantum mechanically at the full MP2/6-31++G(d,p) level of theory. From these data, we have calculated the temperature effects on the cluster thermodynamic functions, and thus the equilibrium Boltzmann distribution for each n. For a selected number of isomers, we have corrected the calculations for basis set superposition error (BSSE) to obtain accurate clustering energies, in excellent agreement with experiment. The computed clusters are overwhelmingly 4-coordinated, as opposed to bulk liquid water, where sodium cations are believed to be mostly 5- or 6-coordinated. To explain this, we suggest the “cluster stability rules”, a set of coordination-number-dependent hydrogen-bond (HB) strengths that can be obtained using a single BSSE correction. Assuming additivity and transferability, these reproduce the relative stability of most of our computed isomers. These rules enable us to elucidate the trends in HB strengths, outlining the major determinants of cluster stability. For n = 4 and 5, we have also performed anharmonic vibrational calculations (VPT2) to compare with available photodissociation infrared spectra of these gas-phase clusters. The comparison suggests that the experiments actually monitor a mixture of predominantly 3-coordinated isomers, which is quite remote from the computed Boltzmann distribution, particularly at low temperatures. Surprisingly, for these experiments, water evaporation pathways can rationalize the non-equilibrium isomer distribution. The equilibrium isomer distribution is, in turn, rationalized by the entropy of internal rotations of “dangling” water molecules. closer to 5 than to 6.15,18−21 The consensus, then, is that the liquid-phase CN is between 5 and 6. Gas-phase clusters of an ion are thought to mimic its local solvation. However, the CN of (small) aqueous sodium ion clusters, Na+(H2O)n, is closer to 4 than to 6 (or 5), and hence “the first solvation shell of Na+ in the gas phase is more tetrahedral in nature, unlike the octahedral configuration adopted in the bulk phase”.22 For a while, this was the subject of intense debate in the theoretical chemistry community.23−30 Computational chemistry provides information on the arrangement of solvent molecules around metal ions, their binding energy, and the vibrational frequencies of the complexes that are often useful for interpreting experiments. The interplay between theory and experiment is crucial for obtaining a molecular-level understanding of ion solvation, which will help in elucidating the correct mechanisms that control ion movement in biology. The determination of the thermodynamic potentials and the assignment of the infrared (IR) spectra of hydrated ions by means of ab initio calculations are subject to the identification of all the low-energy isomers of the cluster at a given temperature (T). We designate these

1. INTRODUCTION Metal ions are ubiquitous in chemical and biochemical processes. The sodium cation, in particular, plays an important role in diverse physiological processes. It is an essential nutrient that regulates body fluids, blood volume, blood pressure, osmotic equilibrium, pH, transmission of nerve impulses, heart activity, and certain metabolic functions in animal and human bodies.1 Its hydration is important for understanding the underlying molecular mechanisms of such processes; for example, the relative selectivity of neural ion channels to sodium vs potassium ions could arise from the need to partially desolvate the ion upon channel entry and replace some of the hydration water by coordination to the selectivity filter moieties.2,3 The coordination number (CN), i.e., the number of solvent molecules in the first solvation shell, is thus a physically important attribute for solvated ions. Liquid-phase experiments4−7 and Monte Carlo (MC) or Molecular Dynamics (MD) simulations8−15 indicated that the CN of sodium cations in liquid water is around 6. However, other X-ray and neutron diffraction experiments8,16−18 on concentrated salt solutions suggest CN = 4 or 5. Recent ab initio MD studies of Na+ in bulk water using density functional theory (DFT) yielded CNs © 2016 American Chemical Society

Received: January 12, 2016 Published: February 25, 2016 1656

DOI: 10.1021/acs.jctc.6b00038 J. Chem. Theory Comput. 2016, 12, 1656−1673

Article

Journal of Chemical Theory and Computation isomers by pqr, where p, q, and r represent the numbers of water molecules in the first, second, and third solvation shells, respectively. Evidently, p is the CN and, in the absence of fourth-shell ligands, p + q + r = n. The n = 6 cluster, in particular, has been studied extensively, with the hope of obtaining better understanding of its hydration in liquid water. Historically, the structure of its lowest energy isomer, even at 0 K, was the subject of much controversy.11,24−30 Early ab initio calculations23,26,27 suggested that the hexa-aqua Na+ cation is an octahedral 600 complex, and this was supported by some MD simulation.31 Perez et al.24 reported MD results in which 600 may not be octahedral but rather built of two parallel staggered equilateral triangles, later observed to have S6 symmetry in restricted Hartree−Fock ab initio calculations.32,33 Lybrand and Kollman25 first claimed that a 420 isomer, with two second-shell water molecules accepting two hydrogen bonds (HBs) each, is the global minimum energy structure of their many-body potential energy functions for Na+(H2O)6. Room-temperature MD (and MC) simulations gave apparently conflicting results. Cieplak et al.28 reported simulations favoring the 510 conformer, while later work29 confirmed the 420 structure. Perera and Berkowitz30 found the latter to be dominant for the SPCE/POL water model, but with the TIP4P water model they obtained mostly 600, with occasional restructuring to give the 510 isomer. These results suggest that the three above-mentioned isomers may be close in energy. In such circumstances, empirical potential functions employed in MD simulations may not be sufficiently accurate to properly determine the energetic isomer order. To resolve this, Kim and co-workers34,35 (extending earlier work by Bauschlicher et al.36) employed the ab initio MP2 method with the [6s5p2d] basis set for Na+, and the DZd and TZ2P basis sets for water molecules, in exploring the various conformers of Na+(H2O)n=4−6 clusters. They have indeed found small energy differences between the 420, 510, and 600 conformers, such that the actual preference is governed by zero-point energy (ZPE) and finite-temperature effects. Thus, while at 0 K 420 is the most stable structure, at 298 K 420 and 510 have comparable structural stabilities. Subsequent ab initio quantum chemistry calculations consistently implicated 420 (D2d symmetry) as the lowest energy n = 6 isomer at 0 K.37−41 In particular, Neela et al.41 performed an exhaustive screening of a large number of Na+(H2O)n=1−20 isomers at the MP2/cc-pVTZ//B3LYP/ 6-31G(d,p) level of theory. They found that, for complexes up to n = 12, CN = 4 is preferred, while CN = 5 and 6 are favored for larger complexes containing up to 20 water molecules. Thus, the transition to liquid-phase coordination occurs at rather large cluster sizes. As noted by Kim et al.,35 the isomer distribution can change markedly with T, due primarily to an entropic effect. The major contributors to the entropy differences between isomers are the many soft vibrational modes that large clusters possess. Once the vibrational frequencies are obtained from an ab initio computation, the temperature-dependent entropy term can be assessed in the so-called harmonic superposition approximation.42 This has been done routinely for, e.g., water clusters43 and charged water clusters,44 and most extensively for protonated water45−50 and methanol clusters.51,52 For the hydrated sodium cation clusters, entropy-driven isomers were investigated (as a function of T) by Miller and Lisy,53 but only for n = 4.

Aside from the usual rigid-rotor and harmonic-oscillator assumptions used to evaluate the rotational and vibrational partition functions, one utilizes the Boltzmann distribution, which is valid only under equilibrium conditions. The question, then, is whether and to what extent gas-phase clusters are equilibrated? The question arises because, in dilute gases at low pressures, there are few collisions that could equilibrate the isomeric distribution, particularly at low T’s. IR measurements on gas-phase hydrated sodium ion clusters that could address this issue were conducted by Lisy and co-workers.53−55 In their infrared Ar-predissociation (Ar-IRPD) experiments, the clusters achieved low T’s (ca. 100 K for n = 4), mainly by successive events of Ar evaporation from the cluster.54−56 In the infrared multiphoton dissociation (IRMPD) setup, water molecules evaporate and the clusters maintain a higher internal T (ca. 300 K for n = 4).53 Thus, these two methods sample the behavior at both low and high T’s, though only in the high-frequency OH stretch regime (3300− 3900 cm−1). Surprisingly, perhaps, one finds that, at low T’s, the most stable isomer dominates only for n = 2 and 3. For n = 4, the experimental spectrum agrees (semi-quantitatively) with the spectrum computed for the 310 isomer, rather than 400 (the isomer of global free energy minimum at all T’s).54 For n = 5, it was concluded that the 311 isomer dominates, in addition to the 410 global minimum.55 Thus, at low T the isomer distribution is not at equilibrium, with enhanced probability for isomers that are more linear and less branched than the global minimum one. Interestingly, the same conclusion was reached by a recent analysis of the IRPD spectra of protonated water clusters, H+(H2O)n, for n = 4 and 5.57,58 Rather (or in addition to) the global minima branched isomers with an H3O+ core, this analysis implicated the involvement of linear isomers with H5O2+ cores. The present ab initio study of the temperature-dependent isomer distribution in Na+(H2O)n clusters extends the work of Kim and co-workers35 and Lisy and co-workers53−55 in several directions: • Studying clusters up to the octamer, n = 8, rather than the customary range of n = 1−6. • Finding some isomers that were not reported before contributing to the isomer distribution, even for n ≤ 6. • Reporting isomer distributions as a function of T for n > 4; data not previously presented. • Comparing IRPD and IRMPD spectra for n = 4 and 5 with anharmonic vibrational spectra computed using second-order vibrational perturbation theory (VPT2).59 The last enhancement allows for a more quantitative comparison of experimental and computed IR spectra, which is necessary for assigning many peaks within a limited frequency range. Last, but not least, the availability of a large body of data on relative isomer stabilities generated herein allows us to embark on a rather detailed investigation concerning the determinants of isomer stability. This, evidently, includes HB energies of water ligands in the first and second hydration shells of the sodium cation, which depend on the CN and the type of the HB. The problem is that binding energies are prone to basis set superposition errors (BSSEs),60 which we note are quite substantial in the first solvation shell. Computing all the isomers with, e.g., the counterpoise (CP) correction of Boys and Bernardi61 is a huge task when so many isomers are 1657


Article

Journal of Chemical Theory and Computation involved. We find that extensive BSSE corrections can be circumvented by having just a single HB energy that is corrected for BSSE. We summarize these findings as “stability rules” that allow, for the first time, to deduce the relative isomer energies and outline the factors determining them.

Under the ideal gas assumption, the enthalpy change is the same as the internal energy change because, at a given T, ΔRT = 0, where R is the ideal gas constant. Hence: ΔG(T ) = ΔU (T ) − T ΔS(T )

(1)

where ΔU and ΔS are the internal energy and entropy differences. The electronic and translational partition functions do not contribute to either function, because all isomers (for a given n) have the same spin state and the same mass. Additionally, in the classical (equipartition) limit, the rotational energy of a non-linear molecule is 3RT/2. As this is identical for all isomers, it does not contribute to ΔU. Subsequently:

2. COMPUTATIONAL DETAILS 2.1. Computational Level of Theory. Geometry optimizations and calculations of energies and harmonic frequencies were carried out, using the Gaussian09 suite of codes,62 at the full MP2/6-31++G(d,p) level of theory for most of the low-energy isomers of Na+(H2O)n=1−8 clusters in the gas phase. For a given n, the initial structures for optimization were obtained systematically by considering all the possible ways of attaching a water molecule to any of the n − 1 isomers. In these optimizations, energies were converged to within 10−7 hartree, while the root-mean-square force and its largest component, the root-mean-square of the displacement and its largest component (in atomic units), were under 0.00030, 0.00045, 0.00120, and 0.00180, respectively. The coordinates of all optimized structures are collected in Table S1 of the Supporting Information (SI). The symmetry of the isomers was determined by visualization. It is approximate because the more stringent cutoff applied in the Gaussian code often results in no symmetry elements (C1 symmetry). Hessian matrices were also calculated in order to verify that the stationary points on the potential energy surfaces are minima (and not transition states) and to calculate harmonic vibrational frequencies and intensities. To compare with the available experimental results, anharmonic frequencies and intensities were calculated for the isomers of Na+(H2O)n=4,5 using VPT2,59 which is available in Gaussian09, Rev. D.62 Results of these calculations for n = 4 and 5 are collected in Tables S2 and S3 of the SI, respectively. BSSE60 can be quite substantial for water attachment to a positively charged cluster. The BSSE arises because, as two monomers approach one another, they can utilize each other’s basis set to lower the dimer’s electronic energy, thus overestimating the attractive interaction in the dimer (namely, the HB strength). However, when isomers with the same number, n, of water molecules are compared, BSSE is not expected to be an issue. Therefore, throughout most of this work we do not correct for BSSE, except toward the end when we also report sequential water binding energies in comparison with experiment.63 There we utilize the CP correction as suggested by Boys and Bernardi.61 We use one of the CP-corrected HB energies to calibrate all other HB energies in the “cluster stability rules” described below. 2.2. Calculating Thermodynamic Parameters. Thermodynamic functions at 298.15 K and 1 atm for all isomers were calculated in Gaussian and collected in Table S1 of the SI. Temperature effects on cluster thermodynamic data were evaluated by running a corrected version of our FORTRAN Tempo code.52 This code uses the same thermodynamic formalism (below) as the Gaussian code but allows for efficient calculations of the cluster thermodynamics as a function of T (rather than at just one specified T as in the Gaussian output). The free energy difference, ΔG, between isomers of a given cluster (same n) was evaluated from the rotational and vibrational partition functions (qrot and qvib, respectively) in the rigid-rotor harmonic-oscillator approximation.52 All computations have been done considering 1 mol of our molecular system as an ideal gas described in the canonical ensemble.

ΔU (T ) = ΔEelec + ΔUvib(T )

(2)

ΔS(T ) = ΔSrot + ΔSvib(T )

(3)

Eelec is the electronic energy of the molecular system (conventionally measured relative to the dissociation limit of non-interacting neutral atoms). It is retrieved from the Gaussian output file and does not depend on T. The subscripts “rot” and “vib” denote rotational and vibrational entities. Because all the clusters are non-linear (at least when hydrogens are included) and the T is not extremely close to 0 K, their rotational partition function can be evaluated in the classical (high-T) limit: π qrot(T ) = σ

⎛ 8π 2k T ⎞3/2 B ⎟ ⎜ I1I2I3 2 ⎠ ⎝ h

(4)

From it, one obtains the rotational entropy change: ΔSrot = R Δ ln qrot(T )

(5)

I1, I2, and I3 are the principal moments of inertia for a nonlinear moiety, obtained from the Gaussian output file. More precisely, Gaussian outputs qrot(T0), at a preselected temperature T0, are scaled by (T/T0)3/2 to obtain qrot(T). For obtaining ΔSrot this scaling is not really necessary, because T3/2 cancels when one takes the ratio of the rotational partition functions of two isomers, leading to a temperature-independent ΔSrot. σ is the symmetry factor (the order of the rotational subgroup of the cluster). Also, kB is Boltzmann’s constant, and h is the Planck constant. For N atoms (N = 3n + 1), there are 3N − 6 vibrational modes. Assuming that no imaginary frequencies are obtained in the calculations, we have 3N − 6

qvib(T ) =

∏ i=1

θi =

e−θi /2T 1 − e−θi / T

(6a)

hνi kB

(6b) 3N − 6

Svib(T ) = R

∑ i=1

3N − 6

Uvib(T ) = R

∑ i=1

⎡ ⎤ θi − ln(1 − e−θi / T )⎥ ⎢⎣ θi / T ⎦ e −1

(7)

⎛1 1 θi⎜ + θ / T i ⎝2 e −

(8)

⎞ ⎟ 1⎠

The (harmonic) vibrational frequencies, νi, with i = 1, ..., 3N − 6, are obtained directly from the Gaussian output. In this formalism, hindered internal rotations are treated as soft vibrational modes. 1658


Article

Journal of Chemical Theory and Computation Note that the ZPE, R∑i θi/2, is included in the last equation for the vibrational energy. It is customary to add it to the electronic energy to obtain the “ZPE-corrected” electronic energy, ′ ≡ Eelec + R ∑ θi /2 Eelec i

water dipole moment vector (the HOH bisector) lies coaxially with the Na+O axis. Assuming that the water molecule retains its HOH angle of 104.5°, the two Na+OH angles are expected to be (360 − 104.5)/2 = 127.8°. In our MP2 calculations, we obtained 104.7° for the HOH angle and 127.7° for the two Na+OH angles, in close agreement with this expectation. The small differences could be explained by the Valence Shell Electron Pair Repulsion (VSEPR) model.64 According to it, the reason that the HOH angle in gas-phase water is smaller than the tetrahedral angle of 109.5° (sp3 hybridization) is repulsion with the oxygen lone-pair (LP) electron density. Now that some of this electron density is transferred to the Na+ cation, there is less repulsion. The water angle is then slightly larger than 104.5°, and consequently the Na+OH angles are slightly smaller than 127.8°. The same coaxial water orientation is seen for larger clusters (below), particularly for first-shell water ligands that are not distorted by further hydrogen-bonding. Interestingly, it also characterizes cation solvation in liquid water, according to fs-IR experiments.65 The coaxial water orientation is not self-evident, because water is often envisioned as having two LP charge densities concentrated at the tetrahedral positions.66 Had the sodium cation been bound to just one of the LPs, the Na+OH angle would have been close to the tetrahedral angle of 109.5°. But in reality there are “no rabbit ears on water”.67 The water LP electron density is distributed between the two tetrahedral sites, along the so-called “negativity track”, and thus there is nonnegligible electron density also at the “trigonal site” (along the HOH bisector) to which the cation binds.68 This allows the sodium cation to interact favorably with the trigonal site electron density, and equally with the two tetrahedral LP sites. 3.2. The Dimer. For n = 2, there is also a single stable isomer for all T’s35−39,41,54,55,69,70 (see Figure 1b). The two water molecules are coordinated with the sodium cation (200 arrangement), and not hydrogen-bonded to each other (unobserved 110 arrangement). This has been verified by Ar-IRPD experiments, which show only the symmetric and asymmetric OH stretch bands of free OH groups, and no trace for a hydrogen-bonded OH at lower frequencies.54 It was explained by the fact that the Na+···(OH2)2 interactions in the first solvation shell are larger than the Na+(OH2)···(OH2) interactions in the second solvation shell. The two water oxygen atoms are diametrically opposed along the ONa+O line, and their hydrogen atoms are staggered so that the two (HOH) water planes are perpendicular to each other, resulting in the D2d symmetry. This arrangement can be rationalized by electron repulsion, as emphasized by the VSEPR model, which in turn is a consequence of the Thompson problem of finding the geometry that minimizes the Coulombic energy of n like charges that are free to move on the surface of a sphere (see http://mathpages.com/home/kmath005/ kmath005.htm). For a given Na+O bond length, the diametrically opposing positions of the negatively charged oxygen atoms keeps them farthest away from each other. In addition, the repulsion is manifested in the Na+O distances, which are slightly larger in the dimer (2.26 Å) than in the monomer (2.25 Å). The staggered hydrogen positions keep the LPs of the two oxygen atoms at their maximal possible distance. Indeed, the LPs point toward the central Na+, so that minimizing their repulsion (rather than the hydrogen atom repulsions) must be dominant in determining the angles between the two water

(9)

In this work, all electronic energies are ZPE corrected (even when not explicitly stated). Let us designate by ΔGj(n) the free energy difference (say, from the most stable isomer at the given T) for the jth isomer in a cluster of size n. The relative stability of the kth isomer (probability P(n) k ) under equilibrium conditions will be given by the canonical probability distribution (n)

Pk(n)(T ) =

e−ΔGk ∑j e

(T )/ RT

−ΔG(j n)(T )/ RT

(10)

Here j runs over all possible isomers. The comparison of this expression with the experimental isomer distribution, when measured, could indicate whether the cluster in the gaseous sample is under equilibrium or not.

3. RESULTS We now discuss the 75 isomers of the Na+(H2O)n=1−8 clusters that were found in this work by minimization of the electronic energy. For each n, the investigated isomers are within approximately 10 kcal/mol of the lowest energy isomer. We present in the figures below their structures, symmetries, and zero-point-corrected electronic energies. Proceeding systematically in the order of increasing n and ΔEelec ′ , new HB types are encountered. Their energies are then recorded in the “stability rules” table. Using it, one may predict the relative stabilities of all the other isomers. For n > 4 we present the temperature dependence of the isomer distribution at equilibrium. For n = 4 and 5 we have also computed the anharmonic (VPT2) spectra, for comparison with the IRPD and IRMPD measurements performed in the OH stretch regime. 3.1. The Monomer. The n = 1 structure has been well known since the first reliable (MP2) quantum chemistry calculations for these clusters.36 It has a single stable isomer with C2v symmetry at all T’s (see Figure 1). This means that the

Figure 1. Isomers of the gas-phase Na+(H2O)n=1−4 clusters. Numbers in blue represent the zero-point corrected electronic energy differ′ (in kcal/mol), of the corresponding isomers from the ences, ΔEelec structure of global minimum in electronic energy for given n. The symmetry point group of each structure is indicated in green. 1659


Article

Journal of Chemical Theory and Computation planes. The inward directionality of the water LPs (more precisely, of its “negativity track”) in aqueous Na+ make the Na+ occupy an internal site, whereas in neutral sodium clusters the sodium atoms tend to occupy a surface position.37 3.3. The Trimer. The n = 3 cluster continues this trend, with a triply coordinated Na+ (isomer 300 in Figure 1c) being the most stable at 0 K,35−39,41,54,55,69,70 as well as at elevated T’s. As Table S4 in the SI shows, this isomer has the lowest free energy at all T’s. Like for the dimer, the Ar-IRPD spectrum for the trimer shows only free OH absorptions, and none for water−water HBs.54 Thus, the first-shell water−water repulsions are still not large enough to push the third water molecule into the second solvation shell. These repulsions are nevertheless manifest in the water oxygen-atom arrangement, which is an equilateral triangle with the sodium cation in its center (D3h symmetry). This arrangement minimizes the Coulomb repulsion of three equal charges that are free to move on the surface of the sphere. Moreover, due to LP−LP repulsions, the OH vectors are tilted from the normal to the trioxygen plane by about 60.0°, so that the H atoms are closer to this plane than the LPs, and the symmetry is reduced from D3h to its rotational subgroup, D3. For n = 3 we have also found two higher energy isomers, 210-1 and 210-2 (see Figure 1d,e), with similar electronic energies, about 5 kcal/mol higher than the 300 isomer. This agrees with Miller and Lisy,54 who have shown that these two isomers do not contribute to the high-temperature IRMPD spectrum.53 The transition from 300 to 210-2 ( ΔE′elect = 5.4 kcal/mol), involves detaching one water molecule from the inner shell and attaching it in the second shell. Thus, one wonders whether the relative isomer stabilization energies can actually be obtained as the difference between two water attachment energies. Unfortunately, substantial BSSEs are involved in the calculation of the clustering energies. Such errors are not involved in the relative isomer stability energies for fixed n. Although the CP correction can remedy BSSEs, this requires additional calculations that are sometimes debatable.60 However, as the example of cluster 210-2 shows, if we knew one of the water attachment energies, we could determine the others from the relative isomer stabilities without the need to invoke the CP correction. This is the basis for our “cluster stability rules” below. Let us denote by ΔEp→p−1 the ZPE-corrected, electronic energy for detaching a first-shell water molecule (W1) when the CN = p, namely, for breaking one of the p Na+···W1 bonds. We have calculated this energy for p = 3 using the CP-correction for the BSSE, obtaining ΔE3→2 = 16.1 kcal/mol. We use this value below as a basis for obtaining all other HB energies. Consider next a HB, W1→W2, donated from W1 to a secondshell water molecule, W2. If W1 and W2 are not involved in any other HB (a “simple” HB) and CN = p, we denote its ZPEcorrected 0 K energy by ΔEw(p). Now, ΔEw(2) for a doubly coordinated sodium cation can also be calculated from the relative isomer energies along the path 300→210-2→200. This gives

Table 1. Stability Rule Electronic Energies (kcal/mol) for Hydrated Sodium Ion Clusters as a Function of the Na+ Coordination Number, pa ΔEp→p−1b

ΔEw(p)c

ΔEAA(p)d

ΔEaDD(p)e

2 3 4 5 6

16.1g 14.6 10.6 11.4

10.7 9.9 8.9 9.4 ”h

0.6 2.3 2.8 2.3 ”h

9.1 8.8 ”h

ΔEsDD(p)f 9.1 8.6 ”h

a

Entries are justified as we proceed from small to large clusters. Thirdshell HBs are assumed to be independent of p. ΔEAD = 8.4 kcal/mol (energy for removing an AD W2 water molecule in the W4→W1←W2← W3 HB network, and forming the W4→W1←W3 network), ΔEAAD = 10.1 kcal/mol (energy for cleaving the W2→W3 HB when W2 is a second-shell water ligand of AAD type), and ΔEAADD = 9.1 kcal/mol (energy for removing a W2→W3 HB when W2 is a second-shell water ligand of AADD type), as well as the second- and third-order coupled AA HBs, ΔEAA2 = 1.2 kcal/mol and ΔEAA3= 0.5 kcal/mol. Second-shell HBs cease to depend on p beyond p = 5. All energies are in kcal/mol and are ZPE-corrected. bInner-shell water molecule (W1) evaporation energy, when W1 is not further hydrogen-bonded and CN = p. c Second-shell water molecule (W2) evaporation energy, when W1 and W2 are not further hydrogen-bonded and the sodium cation is in the coordination p. dEnergy for cleaving one of the W1→W2 bonds to a double-acceptor W2 molecule. eEnergy for cleaving a second W1→W2 HB, when W1 is of DD type in an asymmetric environment. fEnergy for cleaving a second W1→W2 HB, when W1 is of DD type in a symmetric environment. gCalculated using the CP correction; see Table 6, below. hValue as in the preceding row.

A closer inspection shows that the 210-1 isomer is more stable than 210-2, suggesting that the double-acceptor water molecule (AA) in its second shell is stabilized in comparison to the single-acceptor (A) second-shell water molecule of 210-2, by the formation of the second HB. However, because of the large O−Na−O angle change required for forming the AA bond when CN = 2, this extra stabilization is worth only 0.6 kcal/mol. We denote it by ΔEAA(2). The extra stabilization of the AA bond comes from the fact that it is actually two HBs. When broken sequentially, the first is worth ΔEAA(2) and the second is equal to ΔEw(2) . However, when both exist they are equivalent, the HB energy of each being [ΔEw(2) + ΔEAA(2)]/2. Thus, each one is a weak HB, so that in the IR spectrum the OH stretch for the W1 molecule participating in this HB is only slightly red-shifted from the free OH case. 3.4. The Tetramer. 3.4.1. Structure. For n = 4, 14 structures were investigated. After geometry optimization, we retained four different stable isomers (see Figure 1f−i). In agreement with the literature,35−39,41,54,55,69,70 the most stable isomer is 400. Its geometry can also be explained using the VSEPR theory. The solution to Thompson’s problem is now a tetrahedron, so that VSEPR predicts Td symmetry for the four water oxygens, with Na+ in the center. The orientation of the hydrogen atoms again reduces the symmetry to that of a rotational subgroup, S4. We located three higher energy isomers, 310, 220-1, and 220-2 (Figure 1g−i). Of these, only 310 was found in the computations of Miller and Lisy.54 We can use these data to determine the electronic energy for water detachment from the 400 isomer without invoking the CP correction. For this aim, we consider the path 400→220-1→210-2→300. The step 220-1→210-2, in which n changes, involves a previously determined parameter, ΔEw(2) = 10.7 kcal/mol, whereas the

ΔEw (2) = ΔE3 → 2 − ΔE210‐2 = 16.1 − 5.4 = 10.7 kcal/mol

p

(11)

We report this value in Table 1, which we shall complete as we progress to larger n values. 1660


Article

Journal of Chemical Theory and Computation other steps involve relative isomer energies from Figure 1. Thus, we obtain ΔE4 → 3 = ΔE220‐1 + ΔEw (2) − ΔE210‐2 = 9.3 + 10.7 − 5.4 = 14.6 kcal/mol

(12)

This number, although smaller than ΔE3→2, is still substantially larger than ΔEw(3) below, which explains why 400 is the global minimum for n = 4. Consider next the relative stabilities of the remaining two isomers. The energy of the transition 400→310 is ΔE4→3 − ΔEw(3) − ΔEAA(3). We determine below that ΔEAA(3) = 2.3 kcal/mol, and therefore ΔEw (3) = ΔE4 → 3 − ΔE310 − ΔEAA (3) = 14.6 − 2.4 − 2.3 = 9.9 kcal/mol

(13)

As expected, ΔEAA(3) > ΔEAA(2) as the AA bond becomes less strained, and ΔEw(3) < ΔEw(2) as the solvation burden is distributed among more water ligands. The 220-2 isomer is 1.6 kcal/mol less stable than 220-1. The latter has two isolated W1→W2 HBs. Upon conversion to 220-2, one of these is broken, while a second HB is formed to the other W1 molecule. This creates a double-donor (DD) water molecule in the first shell. The DD molecule here is in a symmetric environment, so we call it “symmetric DD” (sDD). Its electronic energy will be denoted by ΔEsDD(p), so that ΔE220‐2 − ΔE220‐1 = ΔEw (2) − ΔEsDD(2)

Figure 2. Comparison of experimental IR spectra53−55 for the tetraaqua sodium ion cluster with anharmonic VPT2 calculations. Black line is the infrared photodissociation (IRPD) spectrum of Na+(H2O)4Ar at T ≈ 100 K,54 whereas the red line depicts the infrared multiphoton dissociation (IRMPD) spectrum of Na+(H2O)4 at T ≈ 300 K.53 Anharmonic spectra of the four isomers of Na+(H2O)4 from Figure 1 (no Ar attached), calculated at full MP2/6-31++G** level of theory, are represented by colored stick diagrams (see Table S2 in the SI). A scale factor of 0.9936 has been applied to all calculated frequencies to match the experimental frequency of the bent (AA) HB stretching mode of the only isomer exhibiting this feature (310). Each isometric spectrum was normalized to unity at its highest peak. Experimental data courtesy of James M. Lisy.

(14)

which gives ΔEsDD(2) = 9.1 kcal/mol. Thus, the energy required to cleave the first of the two W1→W2 HBs is 9.1 kcal/ mol, while the HB energy of the remaining W1→W2 HB is 10.7 kcal/mol. Consequently, a DD water molecule present in a cluster acts as a destabilizer, whereas an AA water molecule acts as a stabilizer. 3.4.2. Spectroscopy. It is interesting to compare the relative isomer stability with spectral evidence for isomers contributing to the low- and high-temperature IR spectrum. Three OH bands were observed in the Na+(H2O)4Ar IRPD spectrum (low T) of Miller and Lisy:54 two non-hydrogen-bonded features, 3715 and 3649 cm−1, assigned to the asymmetric and symmetric stretch, respectively (of the free OH bonds), and a prominent hydrogen-bonded feature at 3547 cm−1. Based on harmonic normal-modes, the latter was assigned to the “bent HB” (i.e., AA) in the 310 isomer. This band is the least redshifted from all the hydrogen-bonded OH stretches, because the two AA bonds are the weakest, each with an average HB energy of [ΔEw(3) + ΔEAA(3)]/2. However, the harmonic calculations are not sufficiently accurate for determining whether the global minimum isomer (400) is present. A more definitive test is provided by the anharmonic (VPT2) IR spectra for the four n = 4 isomers (Table S2 in the SI), depicted in Figure 2 as a stick diagram. In contrast to the large blue-shift of the harmonic frequencies requiring a scaling factor of 0.9604,54 the anharmonic spectra are only mildly blue-shifted, requiring a scaling factor (0.9936) that is much closer to unity. This correction factor could be attributed to the absence of nuclear quantum effects in our calculations. As seen in Figure 2, the VPT2 calculations clearly support the conclusion of Miller and Lisy that the dominant conformer in the low-temperature IRPD spectrum is 310. After scaling, there is still a residual blue-shift of the calculated 3715 cm−1

band and a slight red-shift of the 3649 cm−1 band (blue bars). We ascribe these to the effect of the attached Ar atom, which was not included in our calculations. It red-shifts the asymmetric stretch by 17 cm−1 while blue-shifting the symmetric stretch by 5 cm−1 (compare Table 2 in ref 54 with Table 2 of ref 53). Adding these shifts will bring our calculations to perfect agreement with experiment. Moreover, even the two doublets (somewhat smeared when Miller and Lisy smoothed their data) in the 3715 and 3547 cm−1 bands are reproduced computationally, which leaves little doubt that the IRPD spectrum is indeed that of the 310 isomer. In contrast, the antisymmetric stretching mode of the 400 isomer is red-shifted by 50 cm−1 from the experimental peak. It will not agree with the IRPD 3715 cm−1 band even if we redshift it by the extra 17 cm−1, indicating that 400 does not participate significantly in the Na+(H2O)4 population at low T. The 220 isomers are even less likely participants because they have peaks around 3400 and 3500 cm−1 that do not appear in the IRPD spectrum (Figure 2). Thus, the observed IRPD spectrum can be explained by a nearly pure 310 contribution. We note that the computed hydrogen-bonded OH stretching modes of the 220 isomers are red-shifted from those of the 310 isomer, indicating that these HBs are indeed stronger than the AA bonds. In particular, the band for the 220-1 isomer is redshifted from that of the 220-2 isomer because the W1→W2 bond is strongest when W1 engages in a single HB. The water− water HBs in the 220-2 isomer are weaker due to the destabilizing sDD interaction. The high-temperature IRMPD spectrum (red line) does not have a 3547 cm−1 band, which rules out 310 participation. The experimental asymmetric stretch band now widens to 1661


Article

Journal of Chemical Theory and Computation

To obtain the water binding energy in the second shell for CN = 4, we use another pathway starting with the same isomer, 500→400→410, hence

encompass the 400 frequency, which is thus dominant in the population of Na+(H2O)4. There is also a possibility that a small fraction of the 220-2 isomer also exists in the molecular beam because a small band near 3450 cm−1 does show up in the IRMPD spectrum. The fact that a non-global-minimum energy isomer (310) is dominant at low T can be due to its trapping during cluster formation. Assume that clusters are formed by sequential water loss until their internal energy becomes too small to promote further water evaporation. If 310 or 400 loses one water molecule, it becomes 300, the only isomer observed for the trimer. In contrast, if any of the 220 isomers loses a water molecule, they convert to 210, which is not observed for the trimer. Hence, the 220 isomers cannot exist in the low-T tetramer. This simple argument would, however, be invalid if the clusters cool exclusively by evaporating Ar ligands, as assumed by Miller and Lisy.54 In contrast to the situation at low T, the high-T population appears to be closer to equilibrium, because the predicted dominant isomer (400) is indeed the one observed experimentally. 3.5. The Pentamer. 3.5.1. Structure. For n = 5, 17 structures were investigated, and 7 have been retained as different stable structures (see Figure 3). Among these stable

0 ≡ ΔE410 = ΔE500 + ΔE5 → 4 − ΔEw (4) − ΔEAA (4) (16)

We now use ΔEAA(4) = 2.8 kcal/mol, which will be determined below from isomers of the hexamer and heptamer, to obtain ΔEw(4) = 1.1 + 10.6 − 2.8 = 8.9 kcal/mol. Thus, the strength of second-shell HBs also decreases with p, but less rapidly than for the first solvation shell. We also use the isomer energies in Figure 3 to complete three additional HB strengths for p = 3. The difference between the 320-3 and 320-1 isomers involves the cleavage of one AA bond. Thus, we get ΔEAA (3) = ΔE320‐3 − ΔE320‐1 = 2.3 kcal/mol

(17)

The difference between the 311 and 320-1 isomers is in moving a water molecule from the second to the third shell, converting an AA water molecule to an AAD type. Its electronic energy is given by ΔEAAD = ΔE320‐1 − ΔE311 + ΔEw (3) = 4.0 − 3.8 + 9.9 = 10.1 kcal/mol

(18)

Thus, an AA water molecule is a slightly better HB donor than the A type. We assume that third-shell HBs are (approximately) p-independent. Another new HB type is encountered in the 320-2 isomer, which involves a (first-shell) DD water molecule in an asymmetric environment. Using the path 320-1→320-2, we find that ΔEaDD(3) = ΔE320‐1 − ΔE320‐2 + ΔEw (3) = 4.0 − 4.8 + 9.9 = 9.1 kcal/mol

(19)

Thus, the aDD water molecule is an inferior HB donor compared to the simple A type. A second donor bond (symmetric or not) for a given water molecule is always weaker than the first such bond. The energy terms determined thus far are collected in Table 1. As a check, we use them to calculate the relative energies for all the n = 5 clusters in Table 2. Note that we have determined five parameters from the n = 5 isomer energies: ΔE5→4, ΔEw(4), ΔEAA(3), ΔEAAD, and ΔEaDD(3), using the five equations in this section. Therefore, the test is in determining the stability of the sixth isomer, 311, which can be calculated along the path

Figure 3. Isomers of the gas-phase Na+(H2O)5 cluster (see legend of Figure 1).

conformers, at 0 K, 410 is the global minimum. This result confirms earlier ab initio calculations by several researchers.35,37,39,41,53−55 At 0 K, the 5-coordinated structure (500) is less favorable than the 4-coordinated 410 isomer by 1.1 kcal/mol. We can estimate the evaporation energy of a water molecule from 500 to give the 400 isomer, again without invoking the CP correction. For that end, we locate a pathway involving only steps of known electronic energy differences: 500→230→220-2→400. Summing these energy differences, we obtain

Table 2. Stability Rule Model Electronic Energies for Penta-aqua Sodium Ion Clusters Relative to the 410 Isomer, Compared with Our MP2 Calculations ΔEpqr (kcal/mol) isomer pqr 500 311 320-1 320-2 320-3 230

ΔE5 → 4 = ΔE230 − ΔE500 + ΔEw (2) − ΔE220‐2 = 11.9 − 1.1 + 10.7 − 10.9 = 10.6 kcal/mol (15)

Thus, there is a dramatic decrease in the binding energy for the fifth water ligand in the first shell, and this contributes to the propensity of the sodium cation to be 4-coordinated in small aqueous clusters.

equation eq eq eq eq eq eq

16 20 18 19 17 15

a

model

MP2

1.1 4.0 4.1 4.9 6.4 11.9

1.1 3.8 4.0 4.8 6.3 11.9

a

Only model (and not MP2) energies are used in these equations.

1662


Article


present. According to the thermodynamic stability curves of Figure 4, 410 is the dominant isomer at all T. Thus, the cluster distribution in the high-T experiments is closer to being equilibrated, whereas at low T the system is clearly out of equilibrium. If these assignments are correct, and clusters are formed by water evaporation (in addition to the Ar evaporation process54), one can identify several water loss pathways that contribute to the buildup of the observed isomer population. The first is, evidently,

410→400→300→310→311 using energy terms from Table 1 as follows: ΔE311 = ΔEAA (4) + ΔEw (4) + ΔE4 → 3 − ΔEw (3) − ΔEAA (3) − ΔEAAD

(20)

This gives 4.0 kcal/mol, which is within 0.2 kcal/mol from the calculated MP2 value. The relative isomer population as a function of T is shown in Figure 4. At higher T’s, the global minimum found at 0 K still dominates the population, yet other conformers rise in significance, particularly 500 and 320-3 (the latter has not been located in previous work). It is thus interesting to confront this equilibrium isomer distribution with the IRPD and IRMPD results, by performing anharmonic VPT2 calculations for the various n = 5 isomers.

(21)

311 → 310 → 300

This pathway can be important at low T’s, as it connects the experimentally dominant low-T isomers for n = 5, 4, and 3. Each step involves cleavage of a water−water HB, which is less costly than evaporating an inner-shell water ligand. Evidently, this mechanism does not occur at high-T, because the lowfrequency bands of the 311 and 310 isomers are not observed in the IRMPD spectrum. In addition, one might have (22)

320‐2 → 310

This water loss process is prompted by the relative instability of the DD HB. At higher T’s one may have (23)

410 → 400 → 300

This channel of water loss connects the most stable n = 5, 4, and 3 isomers. It requires high T’s for cleaving the more energetic inner-shell HBs. 3.6. The Hexamer. For n = 6 we have built 35 structures. From these, 14 distinguishable stable isomers were found (see Figure 6). The lowest electronic energy isomer is 420-1, as noted by many theoretical calculations.25,28−30,35,37,39,41 Binding of the sixth water molecule to the C2 410 isomer results in a more symmetric water hydrogen arrangement, increasing the symmetry to D2d. The 600 isomer has an approximate D2h symmetry if the water molecules bind through their trigonal sites (as above), because then the HOH planes of every pair of adjacent water molecules (those having O−Na−O angle of 90°) are perpendicular to each other in order to minimize LP repulsions. This would result in pairs of diametrically connected water molecules lying in the same plane. The optimized geometry is slightly distorted, with the pairs of hydrogen atoms of the diametrically connected water molecules lying 13.8° above and below the plane. This leads to C1 symmetry, contrary to the S6 symmetry reported by Kim et al.35 This difference could be explained by the fact that for geometry optimizations, Kim et al.35 used the Hartree−Fock (HF) level of theory. This does not include electron correlations and therefore may lead to a different optimized geometry as compared to a more correlated method such as MP2. We can now obtain the water detachment energy from the 600 isomer by following the 600→420-1→410→500 path. This is calculated as

Figure 4. Temperature dependence of the relative stability of different conformers of the Na+(H2O)5 cluster, calculated from eq 10 using the partition functions described in section 2.2.

3.5.2. Spectroscopy. An assignment of the IRPD peaks for Na+(H2O)5Ar (at ca. 75 K) has recently been suggested by Ke et al. based on harmonic normal modes.55 The two dominant isomers identified in the low-T spectrum are 311 and 410. The latter is the global minimum, but it does not contribute to any of the bands at 3360, 3430, and 3465 cm−1. Thus, just like for the tetramer, at least one additional isomer must be present in the low-T sample. Some of our VPT2 anharmonic spectra (Table S3 in the SI) are presented in Figure 5 (without any frequency scaling). If we assume that the asymmetric free OH band gets red-shifted by the Ar tag (that we have not modeled here), then indeed this spectrum confirms the conclusions of Ke et al.,55 except that we find two more isomers that may contribute, 320-2 and 230, which have not been previously identified. Of interest is the 500 isomer, the next in line after the 410 global minimum. While it has not been considered by Ke et al.,55 this omission is likely not crucial for assigning the IR spectrum, because its asymmetric OH stretching band is too blue-shifted in comparison with experiment. At higher T’s, 311 no longer contributes, as its characteristic 3360 cm−1 peak is not manifested in the IRMPD spectrum (red line). The other two isomers, 410 and 320-2, may still be

ΔE6 → 5 = −ΔE600 + ΔEAA (4) + ΔEw (4) + ΔE500 = 11.4 kcal/mol

(24)

This follows by taking for ΔEAA(4) the average of ΔE420‑4 − ΔE420‑1 = 2.7 kcal/mol and ΔE421‑2 − ΔE421‑1 = 2.9 kcal/mol. That is, ΔEAA (4) = 1663

1 (ΔE420‐4 + ΔE421‐2) = 2.8 kcal/mol 2

(25)


Article


Figure 5. Comparison of experimental IR spectra53−55 for the penta-aqua sodium ion cluster with anharmonic VPT2 calculations. Black line is the infrared photodissociation (IRPD) spectrum54 of Na+(H2O)5Ar at T ≈ 75 K, whereas the red line depicts the infrared multiphoton dissociation (IRMPD) spectrum53 of Na+(H2O)5 at T ≈ 300 K. Anharmonic spectra of the possibly significant isomers of Na+(H2O)5 from Figure 3 (no Ar attached), calculated at full MP2/6-31++G** level of theory, are represented by colored stick diagrams (see Table S3 in the SI). No scaling factor was applied to the computed frequencies. The intensities of each isometric spectrum were normalized to unity at the highest peak. Experimental data courtesy of James M. Lisy.

From the hexamer data we find three more p = 4 energies. First,

Again, we assume that third-shell HBs are not p-dependent. To obtain the water binding energy in the second shell for CN = 5, we use the pathway 600→500→510. Hence,

ΔEaDD(4) = ΔEAA (4) + ΔEw (4) − ΔE420‐5 = 2.8 + 8.9 − 2.9 = 8.8 kcal/mol

(26)

ΔEAA (5) + ΔEw (5) = ΔE600 − ΔE510 + ΔE6 → 5

Thus, for p = 4, the energetic difference between ΔEaDD and ΔEw is already negligible. Next, in the isomer 420-2 we find an interesting W-shape formed of four HBs, involving two secondshell AA water molecules that share a common DD HB donor (in the center of the W). Consider forming this shape by sequential addition of HBs, the final step being 420-5→420-2. This step involves forming a second AA HB, denoted AA2, that is weaker than the first AA bond due to this coupling: ΔE AA 2 = ΔE420‐5 − ΔE420‐2 = 1.2 kcal/mol

= 1.4 − 1.1 + 11.4 = 11.7 kcal/mol

Now, from the heptamer data below, we will determine that ΔEAA (5) = 0.5 × ΔE520‐2 = 2.3 kcal/mol

(30)

and hence ΔEw(5) = 9.4 kcal/mol. Finally, in isomer 330-2 there is a sDD HB for p = 3, whose energy is calculated using the path 420-1→600→500→ 320-1→330-2 to give

(27)

Due to scarcity of data, we assume this value is p-independent. In isomer 312 we encounter for the first time a water molecule with four HBs, which we label AADD. Using the path 420-1→600→500→311→312, we then get

ΔEsDD(3) = ΔE600 + ΔE6 → 5 + ΔE320‐1 − ΔE500 − ΔE330‐2 = 1.4 + 11.4 + 4.0 − 1.1 − 7.1 = 8.6 kcal/mol

ΔEAADD = ΔE600 + ΔE6 → 5 + ΔE311 − ΔE500 − ΔE312 = 1.4 + 11.4 + 3.8 − 1.1 − 6.4 = 9.1 kcal/mol

(29)

(31)

Thus, for the same p, the sDD HB is slightly weaker than the aDD HB.

(28) 1664


Article


Table 3. Stability Rule Electronic Energies (in kcal/mol) for Hexa-aqua Sodium Ion Clusters Relative to 420-1, Compared with Our MP2 Calculations ΔEpqr isomer pqr 510 600 411 420-2 420-4 420-5 312 330-1b 420-6c 330-2 330-3d 240e

equationa eq 29 eq 24 ΔEAA(4) + ΔEw(4) − ΔEAAD eq 27 eq 25 eq 26 eq 28 2[ΔEw(4) + ΔEAA(4)] + ΔE4→3 − 2ΔEw(3) − ΔEAA(3) − ΔEaDD(3) 2EAA(4) eq 31 ΔE330‑1 + ΔEw(3) − ΔEaDD(3) 2[ΔEw(4) + ΔEAA(4)] + ΔE4→3 + ΔE3→2 − 2[ΔEw(2) + ΔEsDD(2)]

model

MP2

1.1 1.4 1.6 1.7 2.8 2.9 6.4 6.8

1.1 1.4 1.6 1.7 2.7 2.9 6.4 6.6

5.6 7.2 7.6 14.5

6.9 7.1 7.3 14.5

a

Only model (and not MP2) energies are used in these equations. Calculated using the path 420-1→300→330-1. cCalculated using the path 420-1→420-6. dCalculated using the path 330-1→330-3, because the segment 420-1→330-1 has already been calculated. e Calculated using the path 420-1→200→240. b


We can now test the stability rules of Table 1 in comparison with the MP2 energies in Figure 6. The results are shown in Table 3 for 12 out of the 14 n = 6 isomers (420-1 defines the zero of the energy scale while 420−3 has rather complex HB pattern not covered by the rules). Of these, the MP2 energies of 7 isomers were used in eqs 24−31 for entering new values to Table 1. As a check, these equations are now reused, with model rather than MP2 energies, to recalculate the energies of these isomers. However, the real test for the stability rules is provided by the remaining 4 isomers, whose energies have not been utilized in this process. The detailed equations in the table show how to calculate them using only predetermined terms from Table 1. For three of these isomers the calculated stability rules values are in excellent agreement with the MP2 values (within 0.2 kcal/mol or better). One isomer, 420-6, is an unexplained exception. One can now appreciate that the stability order of the first three isomers for n = 6 (a 4-coordinated, a 5-coordinated, and then a 6-coordinated Na+ cation) is dominated by second shell AA interactions. The 4-coordinated has two such interactions, the 5-coordinated has one, and the 6-coordinated has none. These differences overshadow the differences between the firstshell binding energies,

Figure 7. Temperature dependence of the relative stability of different isomers of the gas-phase Na+(H2O)6 cluster.

population from 0 K and up to 200 K, while the 4-coordinated isomers 411 and 420-5 dominate at higher T’s. These two isomers have not been reported previously in the literature. The 5-coordinated 510 isomer, although only 1.1 kcal/mol above the global minimum isomer, contributes only marginally ( ΔEp + 1 → p

(p = 4 and 5) dictating the observed order. The computed equilibrium isomer distribution as a function of T is shown in Figure 7. The isomer 420-1 dominates the 1665


Article

Journal of Chemical Theory and Computation isomers, respectively. These might be important in explaining the pentamer spectrum as discussed above. Thus, by a water loss mechanism, the dominant n = 6 isomers give rise to the dominant n = 5 isomers. 3.7. The Heptamer. For n = 7, 40 possible structures were built. After geometry optimization, 18 distinguishable isomers were retained, Figure 8. Two global minima were found, 520-1

Table 4. Stability Rule Electronic Energies (in kcal/mol) for Hepta-aqua Sodium Ion Clusters Relative to the 421-1 Isomer, As Compared with Our MP2 Calculations ΔEpqr isomer pqr 520-1

b

430-1 430-2 430-3 610 412c 421-2 421-3 430-4 430-5 430-6 430-7d 520-2 421-4e 340f 331-1g 331-2

equationa

model

MP2

2[ΔEw(4) + ΔEAA(4) − ΔEw(5) − ΔEAA(5)] + ΔEAAD − ΔE5→4 eq 33 ΔEAAD(4) − ΔEaDD(4) eq 35 eq 32 ΔEAA(4) + ΔEw(4) − ΔEAADD ΔEAA(4) ΔEAA(4) + ΔEw(4) − ΔEaDD(4) ΔEAA(4) + ΔEAAD − ΔEaDD(4) ΔE430‑4 ΔE430‑4 ΔE430‑4 + ΔEw(4) − ΔEaDD(4) eq 30 ΔEAAD + ΔEAA(4) + ΔEw(4) − ΔEaDD(4) − ΔEAD ΔEAAD + 2[ΔEw(4) + ΔEAA(4)] + ΔE4→3 − 2[ΔEw(3) + ΔEaDD(3)] − ΔEAA(3) ΔE340 + ΔEw(3) − ΔEAAD ΔE340 + ΔEw(3) − ΔEAD

−0.5

0.0

0.8 1.3 1.7 2.1 2.6 2.8 2.9 4.1 4.1 4.1 4.2 4.6 4.6 7.8

1.2 1.5 1.7 2.1 2.7 2.9 3.2 3.9 4.1 4.2 4.3 4.5 5.0 7.6

7.6 9.2

7.9 9.3

≤0.2

0.0

MAD a

Only model (and not MP2) energies are used in these equations. Calculated using the path 421-1→400→520-1. cCalculated using the path 421-1→412. dPath to 430-4 evaluated above. eCalculated using the path 421-1→421-4. fCalculated using the path 421-1→300→340. g Path to 340 evaluated above. b

which is identical to the n = 5 case. Apparently ΔEw(p) stops changing when p ≥ 5. The transition 430-2→430-1 involves formation of an AA HB (denoted AA3) that connects two V-shapes harboring the two AA bonds in 430-2. Thus, ΔE AA3 = ΔE430‐2 − ΔE430‐1

(33)

This energy difference is tiny (0.3 kcal/mol)much smaller than ΔEAA or ΔEAA2. A similar HB forms in the transition 440-2→440-1 for the octamer (section 3.8, below). There the energy difference is 0.7 kcal/mol. In the sequel we take the average of these two values ΔE AA3 = 0.5 kcal/mol +

Figure 8. Isomers of the gas-phase Na (H2O)7 cluster (see legend of Figure 1).

and assume (due to data scarcity) that it is independent of p. We also find here for the first time the AD water type (in isomer 430-3). To determine ΔEAD, we follow the path 421-1→430-3. Thus,

and 421-1, with C1 and Cs symmetries, respectively. The 520-1 isomer was also located as the global minimum by Neela et al.,41 who have missed the 4-coordinated 421-1 isomer and thus listed 430-1 as the second most stable isomer. 520-1 has the lowest electronic energy also according to the cluster stability rules (Table 4, below). We have found only one 6-coordinated isomer, which we use to deduce the W1→W2 HB energy along the path 610→ 421-1→420-1→600. This gives

ΔEAD = ΔEAAD − ΔE430‐3 = 10.1 − 1.7 = 8.4 kcal/mol (35)

This is the energy gain resulting from inserting a water molecule between two hydrogen-bonded ones. The insertion adds a second-shell donor HB that is seen to be weaker than a first-shell donor HB (aDD or sDD). We assume that it is p independent. Alternately, we can determine from this the AA HB strength between two second-shell water molecules (2.3 kcal/mol), which we do not use in the sequel. We test the stability rules energies (Table 1) against our MP2 calculations in Table 4. The mean absolute deviation

ΔEw (6) = −ΔE610 + ΔEAAD + ΔE600 = −2.1 + 10.1 + 1.4 = 9.4 kcal/mol

(34)

(32) 1666


Article

Journal of Chemical Theory and Computation (MAD) of the 16 entries in the table is below 0.2 kcal/mol, which is indeed an encouraging result, particularly because only four of these entries were used above to generate new parameters for the stability rules. The temperature dependence of the equilibrium isomer distribution for n = 7 is shown in Figure 9. Although the 5-coordinated isomer (520-1) competes strongly with the 4-coordinated isomer 421-1 at low T’s, other tetra-coordinated isomers (412 and 421-3) dominate at higher T’s in such a way that the hepta-aqua sodium cation in the gas phase is exclusively 4-coordinated at T > 200 K. These two isomers have not been reported previously in the literature.

Figure 9. Temperature dependence of the relative stability of different conformers of Na+(H2O)7 (full lines). Dashed lines are the same functions when the rotational partition functions are neglected (i.e., the free energy differences are calculated only from the vibrational partition functions).

We note, again, that the populated isomers of n = 7 are connected with those of n = 6 through the loss of a single water molecule. Thus, when 421-1, which is the dominant isomer for n = 7 at low-to-intermediate T’s, loses a third-shell water molecule, it becomes 420-1, which is the dominant low-T isomer for n = 6. If, however, it loses the second-shell AA water molecule it becomes 411, which dominates the n = 6 distribution at high T. The same isomer is generated also when 412, the dominant n = 7 isomer at high T, loses a third-shell water molecule. The 420-1 isomer is also generated when 520-1 loses an inner-shell water molecule, whereas if it evaporates a second-shell water molecule, the minority species 510 is obtained. 3.8. The Octamer. For n = 8, 33 possible structures were investigated. After geometry optimization, 27 different isomers were identified, Figure 10. Unlike n = 4−7, for which a 4-coordinated isomer was the global minimum, here we find that the hexa-coordinated isomer 620-1 is the global minimum. This contrasts with Neela et al.,41 who find it as the third isomer energy-wise, while the 440-1 isomer, which we find at a higher energy, is their global minimum. This disagreement could occur because they have optimized their structures with DFT (using the B3LYP functional), whereas we have optimized in full MP2. The observation of a CN = 6 isomer as the global minimum occurs here for the first time, and this could be the beginning of an approach toward the liquid-phase behavior. We proceed with the stability rules for the octamer without adding any new parameters to Table 1. We have already noted


that ΔEw(6) = ΔEw(5) = 9.4 kcal/mol, so we can assume that there is no p dependence to the water−water HBs beyond p = 5. In particular, we assume that ΔEAA(6) = ΔEAA(5) = 2.3 kcal/mol. In a similar vein, we assume that ΔEaDD(5) = ΔEaDD(4) and that ΔEsDD(4) ≈ ΔEsDD(3). There are not enough data to determine these values more accurately. Thus, we close Table 1 with five p-dependent parameters (16 total) 1667


Article

Journal of Chemical Theory and Computation and five p-independent ones. With these values, we are now able to derive the stability rule model energies for the octamer, as compared with our MP2 calculations (see Table 5). The results are encouraging since the MAD is just 0.45 kcal/mol, while none of the 26 entries was used for deriving the stability rules parameters in Table 1. The temperature dependence of the equilibrium isomer distribution is shown in Figure 11. From the plot it is evident that the average CN decreases with increasing T. The Table 5. Stability Rule Electronic Energies (in kcal/mol) for Octa-aqua Sodium Ion Clusters Relative to the 620-1 Isomer, As Compared with Our MP2 Calculations ΔEpqr isomer pqr

equationa b

model

MP2

2.1 1.3 2.1

0.7 1.3 1.8

2.3 2.8 3.0 >2.6g

2.4 2.4 2.5 3.1

3.4

3.4

530-1 521c 530-2

same as 530-2 ΔE6→5 − ΔEAAD ΔE6→5 − ΔEaDD(5) − ΔEAA3

620-2d 422-1e 530-3 431-1f

ΔEAA(6) ΔE6→5 + ΔE5→4 − ΔEAAD − ΔEAADD ΔE6→5 − ΔEAD ΔE422‑1 + ΔEAADD − ΔEaDD(4) − ΔEAA3

440-1

ΔE422‑1 + ΔEAADD + ΔEAAD − 2[ΔEaDD(4) + ΔEAA3]

530-4h

ΔE530‑2 + ΔEAA2

3.3

3.5

620-3 440-2i

2ΔEAA(6) ΔE440‑1 + ΔEAA3

4.6 3.9

4.1 4.1

440-3j 440-4k 440-5l 530-5m 422-2n 530-6o 530-7p

ΔE440‑2 + ΔEaDD(4) − ΔEAD ΔE422‑1 + ΔEAADD + ΔEAAD − 2ΔEaDD(4) ΔE440‑4 ΔE530‑3 + ΔEAA(5) ΔE422‑1 + ΔEAA(4) ΔE521 + ΔEAAD + 2ΔEAA(5) − ΔEw(5) ΔE521 + ΔEAAD + ΔEAA(5) + ΔEw(5) − 2ΔEaDD(5) ΔE422‑1 + ΔEAADD + ΔEAA(4) − ΔEaDD(4) ΔE422‑1 + ΔEAA(4) + ΔEw(4) − ΔEaDD(4) ΔE440‑5 + ΔEAA(4) ΔE440‑6 + ΔEaDD(4) − ΔEw(4) ΔE440‑7 + ΔEw(4) − ΔEAD ΔE440‑7 + ΔEaDD(4) − ΔEsDD(4) 2ΔEAA(6) + 2ΔEw(6) + ΔE330‑3 − ΔE600 − ΔEw(3) − ΔEsDD(3) ΔE440‑8 + ΔEAA(4) + ΔEw(4) − ΔEsDD(4)

4.3 4.4 4.4 5.3 5.6 6.6 5.5

4.2 4.4 4.4 4.8 4.9 5.3 5.6

5.9 5.7 7.2 7.1 7.6 7.3 11.1

5.6 6.0 6.6 6.8 7.0 7.2 10.7

10.4

11.0

0.4

0.0

431-2 422-3 440-6q 440-7 431-3 440-8 350r 440-9

MAD

Figure 11. Temperature dependence of the relative stability of different Na+(H2O)8 isomers.

6-coordinated conformer (620-1), the global minimum conformer at 0 K, dominates the population up to 160 K and competes strongly with the 5-coordinated (521) and 4-coordinated (422-1) isomers around 170 K. At room temperature, the population of n = 8 is exclusively dominated by the 4-coordinated isomers. Therefore, to capture the solution-phase hydrated sodium ion structure at high T’s, clusters with (much) more than eight water molecules should be investigated. 3.9. Successive Clustering Energies. To check the reliability of our results presented in the previous sections, we undertook the calculations of the successive clustering electronic energies (ΔEn−1,n), enthalpies (ΔHn−1,n), and free energies (ΔGn−1,n), in the reaction Na +(H 2O)n − 1 + H 2O → Na +(H 2O)n

(36)

Our goal in this section is to compare these results with previous calculations,32,33,37,39−41 with experimental clustering enthalpies,63 and with the outcome from our cluster stability rules. Because the exact isomer distribution in the experiment was not known, previous calculations were conventionally performed for the lowest electronic energy isomers. Specifically, for n = 2−6 these are the 200, 300, 400, 410, and 420 isomers. However, a more meaningful comparison with the room temperature experiments of Džidić and Kebarle63 would involve the most probable isomers at 298 K. According to our calculations, there is no temperature dependence up to n = 5, whereas for n = 6 the most probable room temperature isomer is 411. We have corrected our calculations of these 5 clusters for BSSE, applying the CP correction as suggested by Boys and Bernadi.61 This means that for each n, we have calculated the electronic energy of the n − 1 cluster with an extended basis set that includes “ghost orbitals” for the nth water molecule at its position in the nth cluster.60 Table 6 shows our calculated ΔEn−1,n values for the isomers with the lowest electronic energies, uncorrected and corrected for BSSE, in comparison with the literature. The MADs between the various calculations and our BSSE-corrected results are low. Our results are closest (MAD = 0.5 kcal/mol) to those of Rao et al.,39 who have used the (BSSE-corrected) CCSD(T) method that is usually considered to be the more accurate method. It can be seen that the BSSE is substantial for

a

Only model (and not MP2) energies are used in these equations. Assuming an AA HB not picked up by the plotting software. c Calculated using the path 620-1→600→500→521. dCalculated using the path 620-1→620-2. eCalculated using the path 620-1→400→4221, with p-invariant ΔEw(p) + ΔEAA(p). fThe path 620-1→422-1 was evaluated above. gAssuming there are two AA HBs not picked up by the plotting software. These must be weaker than normal, hence the inequality. hThe path 620-1→530-2 was evaluated above. iThe path 620-1→440-1 was evaluated above. jThe path 620-1→440-2 was evaluated above. kThe path 620-1→422-1 was evaluated above. lSame as 440-4. mThe path 620-1→530-3 was evaluated above. nThe path 620-1→422-1 was evaluated above. oThe path 620-1→521 was evaluated above. pThe path 620-1→521 was evaluated above. qThe path 620-1→440-5 was evaluated above. rCalculated using the path 620-1→600→330-3→350. b

1668


Article

Journal of Chemical Theory and Computation Table 6. Clustering Electronic Energies, −ΔEn−1,n (in kcal/mol), for the Lowest Energy Isomers of Hydrated Sodium Ion Clusters previous worka n

Hashimoto

1 2 3 4 5 6 MADf a

37

Rao

39

shows very low MAD, of about 0.8 and 0.7 kcal/mol for ΔHn−1,n and ΔGn−1,n, respectively. The MD values are also very low, about 0.2 and −0.2 kcal/mol, respectively. Thus, our proposed approach involving the various suggested conformers, and the computational level of theory is adequate to capture efficiently the structures/energies of the hydrated sodium cation clusters, at least around room temperature.

this work 41

Neela

no CP

b

CPc

26.3 24.2 20.0 16.9 16.6 16.6

22.5 20.1 16.8 13.5 12.1 12.0

22.9 19.3 16.4 15.6 12.6 11.0

24.1 22.3 18.5 16.5 13.2 12.9

22.2 19.7 16.1e 13.4 12.9

3.9

0.0

0.9

1.8

0.5

modeld

4. DISCUSSION This paper presents a detailed investigation of structures of hydrated sodium ion, for a large number (75) of isomers of the Na+(H2O)n=1−8 gas-phase clusters. The obtained structures were subject to geometry optimization and energy calculations at full MP2/6-31++G(d,p) level of theory. Their relative stabilities at 0 K were discussed in terms of a set of cluster “stability rules”. In addition, the calculation of (harmonic) normal-mode frequencies for each isomer allowed to compute its thermodynamic functions and, from them, the temperature effects on the distribution of their various conformers at a given cluster size. We proceeded beyond the harmonic approximation, to calculate the anharmonic IR spectra for n = 4 and 5, in comparison with available experimental IR photodissociation spectra. Finally, we have repeated our calculations with CP correction on the lowest 298 K free energy isomer for each n to obtain clustering energies of hydrated sodium ion in comparison with experimental values. The main results of these different lines of analysis can be summarized as follows. 4.1. Cluster Stability Rules. The “cluster stability rules” encompass a set of HB energies that can explain the computed relative stability of most of the isomers for all n. Two underlying principles upon which these rules are based are additivity and transferability. This implies that the energies of the various types of HBs add up to give the isomer electronic energy differences, and that a given type of HB has the same energy in all clusters where it appears. A third principle concerns with the treatment of dissociative bonds and the BSSE. For understanding its implications on the stability rules, it is instructive to differentiate “dissociative energy” (DE) from “non-dissociative energy” (NDE). DE is the energy difference between two particles at infinite separation and at a finite distance (at a local minimum of their interaction potential). NDE is measured between two local minima of a set of N atoms/molecules that do not involve moving a particle to infinity. Thus, the cluster energies that we

16.1 14.6 11.7 11.7 0.6

37

Hashimoto et al. performed calculations at the MP2/6-31+G(d)// HF/6-31+G(d) level of theory with no BSSE corrections. Rao et al.39 used various theory levels, and we present here their most accurate calculations namely, at the CCSD(T)/6-311++G(d,p) level of theory with BSSE corrections (from MP2). Neela et al.41 undertook their calculations at the MP2/cc-pVTZ//B3LYP/6-31G(d,p) level of theory without BSSE corrections. bSum of electronic and ZPE at the fullMP2/6-31++G(d,p) level of theory, see SI. cBSSE-corrected data. d Stability rules, Table 1. eValue adopted for the cluster stability rules, see Table 1. fMAD from the BSSE-corrected CCSD(T) results of Rao et al.39

n ≤ 4, with the MP2 results without the CP correction being on average 2 kcal/mol too high. For n ≥ 5 the added water molecule enters the second solvation shell, and this may be the reason that the BSSE becomes significantly smaller. We can also compare the CP corrected values with those from our stability rules. Up to n = 4 an inner-shell water molecule is evaporated, so that −ΔEn−1,n ≡ ΔEn→n−1. For n = 5 and 6, an AA water is evaporated, at the energetic cost of ΔEw(4) + ΔEAA(4) = 11.7 kcal/mol. In addition, for n = 6 we have calculated the CP-corrected clustering energy for its most stable room-temperature isomer, 411, as 10.8 kcal/mol. This corresponds to cleavage of an AAD HB, whose energy according to Table 1 is 10.1 kcal/mol. The MAD between the stability rules values and Rao’s BSSE-corrected calculations39 is 0.6 kcal/mol, whereas the mean deviation (MD) is only 0.08 kcal/mol. This remarkable agreement confirms the assumptions behind the cluster stability rules. We have also calculated ΔHn−1,n and ΔGn−1,n from the same CP corrected MP2 runs, for comparison with experimental results of Džidić and Kebarle.63 The comparison in Table 7

Table 7. Clustering Enthalpies and Clustering Free Energies of Hydrated Sodium Ion Clusters at Full MP2/6-31++G(d,p) at 298 K (in kcal/mol)a −ΔHn−1,n

−ΔGn−1,n

n

Hashimoto37

Glendening32

Merrill33

Biring40

expt63

this workb

Biring40

expt63

this workb

1 2 3 4 5 6 MAD MD

24.2 22.1 18.0 14.9 14.0 14.0 1.8 1.8

23.5 18.9 16.0 12.4 11.5 10.1 0.7 −0.7

21.8 19.0 16.4 14.8 12.2 11.9 1.0 −0.1

21.6 18.8 15.6 13.0 12.8 12.4 1.1 −0.4

24.0 19.8 15.8 13.8 12.3 10.7 0.0 0.0

23.0 20.8 15.0 14.6 12.8 11.5 0.8 0.2

15.6 12.9 9.6 7.4 0.4 1.9 1.4 −0.9

17.6 13.2 9.3 6.3 3.9 2.9 0.0 0.0

16.9 14.3 8.3 6.7 3.5 2.0 0.7 −0.2

a Hashimoto et al.37 used the MP2/6-31+G(d)//HF/6-31+G(d) level of theory with no BSSE corrections. Glendening et al.32 performed calculations at the MP2/6-31+G*//RHF/6-31+G* level of theory with BSSE corrections. Merril et al.33 undertook calculations at various theory levels, and we retained here their closest results to experimental clustering enthalpies,63 that is, the MP2/aug-cc-pVDZ//EFP/6-31+G* with BSSE corrections. Biring et al.40 performed calculations at the MP2/aug-cc-pVDZ level of theory without BSSE corrections. bFor the most probable isomers at 298 K, namely 100, 200, 300, 400, 410, and 411.

1669


Article


to the tetrahedral geometry of the 400 isomer, in which the angles are just right for binding an additional water molecule in tetrahedral geometry (though the Na+···W distances are shorter than the W···W ones). Note, however, that the total double-HB strength, ΔEw(p) + ΔEAA(p) ≈11.7 kcal/mol, is nearly independent of p. It is only its partitioning between the first and second HB cleavage events that exhibits this p-dependence. These results explain why, for small Na+(H2O)n clusters, CN = 4 is so prevalent. The combined effect of a minimum of ΔEp→p−1 at p = 5 with a maximum of ΔEAA(p) at p = 4 encourages the dissociation of a water molecule from the inner shell of the pentamer and its redeployment as an AA ligand in the second shell of the tetramer, a process that according to Table 1 is downhill by 1.1 kcal/mol. 4.2. Thermal Isomer Distributions. We have calculated the free energies of all the isomers in the rigid-rotor harmonicoscillator approximation. From it we deduced the Boltzmann distribution for the Na+(H2O)n isomers as a function of T, see eq 10. Several new isomers that populate significantly at T ≠ 0 K were included in this analysis. These isomeric distributions show that the coordination number of the hydrated sodium ion is 4 for n = 4−7 at all T’s, in agreement with the cluster stability rules discussed above. For n = 8, CN = 6 at low T (T < 150 K) and between 6 and 4 around 170 K; above this temperature, CN = 4. This could indicate that the liquid-phase behavior of the hydrated sodium ion starts to be mimicked at low T for n = 8. However, to capture the bulk hydrated sodium ion structure at high T, more than 8 water molecules are required. The distributions also show that for n ≤ 5 there is one dominant isomer at all T’s (the global minimum free energy structure). For larger n, the isomer(s) that dominate at low T decrease sharply in population with increasing T while other isomers rise. One wonders what is the origin of this behavior, and which of the many isomers will have enhanced population with increasing T? Clearly this is an entropic effect, so the question is which factors dominate the ΔS values as calculated from the thermodynamic theory (section 2.2). We have already noted that the translational entropy does not contribute to the isomeric entropy difference, since isomers with the same number of water molecules have the same mass. To asses the role of the rotational entropy, we have plotted in Figure 9 the isomer distribution for n = 7, as an example, with (full lines) and without (dashed lines) the contribution from (rigid body) rotation. It is seen that rotation makes some contribution near T = 0, which rapidly diminishes with increasing T. Thus, vibrational entropy dominates the temperature dependence. As T increases, vibrational modes with higher frequencies start contributing, and the relative rotational contribution diminishes. Figure 12 shows the frequency as a function of the vibrational mode number for n = 6 (as an example). The high-frequency intramolecular modes do not contribute to the entropy at room temperature and below. The main contribution comes from the low-frequency “soft” modes, which involve intermolecular motions. Because room temperature corresponds to 200 cm−1, the first dozen of these make the largest contribution. Isomers having lower frequencies than others in this range will have a larger vibrational entropy. As an example, consider the 600 isomer. Its frequencies (cyan line) are higher than the rest, so it has a minimal entropic contribution. Consequently, it does not contribute to the isomer distribution in Figure 7 in spite of having the third lowest electronic energy (Figure 6c).

calculated with Gaussian (listed in the SI) are DEs, but the relative cluster stabilities at fixed n (Figures 1, 3, 6, 8, and 10) are NDEs. DEs are prone to BSSE; NDEs are not. This is the third underlying principle. This principle means that if the absolute energy (relative to dissociated particles) of one isomer is corrected for BSSE (say, by an additive constant α), the absolute energies of all other isomers (relative to the same dissociated state) will also increase by α, while the NDEs will not change. Hence α needs to be calculated only once. We have calculated it for the process 300→200 to obtain ΔE3→2, and this DE set the BSSE correction for all other DEs involved. The HB types in the cluster stability rules (Table 1) can also be divided into these two classes. All of them are DEs except the double acceptor HBs (AA, AA2 and AA3). When one of these is cleaved, the water molecule in question remains bonded by the second A bond, and does not dissociate to infinity. Hence AA energies are NDEs. Consider the equations for relative cluster energies detailed in Tables 2−5. Because they calculate NDEs, there can be any number of NDE terms in any such equation, but the number of positive and negative DE terms must equal (so that α cancels). One can check that this indeed holds true for all these equations. Starting from n = 2 and progressing to larger cluster sizes, we have depositing new HB stability values in Table 1, which presently encompasses 21 different values. These explain the relative energies of nearly all the isomers studied herein up to n = 8. Therefore, it is of interest to discuss the trends in HB strengths revealed in Table 1. Evidently, the strongest interactions are in the first shell, their strengths decreasing as more water molecules are added therein (increasing p) due to water−water repulsion. The surprise, however, is that this trend is not monotonic, with ΔEp→p−1 reaching a minimum at p = 5. It then increases for n = 6, but we have no data for larger p to follow this further. This nonmonotonic behavior has apparently been overlooked in the literature, because the quantity typically calculated there is the evaporation energy (or enthalpy) of a water molecule from the global minimum structure, e.g., Table 6. This structure is n00 up to n = 4. For n = 6, the minimum electronic energy isomer is 420, not 600, and its evaporation energy (420→410) is smaller than ΔE6→5, creating the impression of a continued descending trend. “Simple” HBs in the second shell are, as expected, weaker than those in the first shell; see the ΔEw(p) values in Table 1. Their decrease with p is again non-monotonic, showing a shallow minimum at p = 4. For p ≥ 5, ΔEw(p) no longer changes with p. “Simple” HBs in the third shell, W2→W3, are even weaker (the AD HB), and are sufficiently far from the positive center to be considered as p-independent. If W2 is a double acceptor, its donor capabilities strengthen, so that the AAD and AADD HBs are actually stronger than the AD HB. In water clusters, AD water−water HBs are stabilized due to cooperativity, leading to ring structures in which all water molecules are of the AD type.71 In the second solvation shell of Na+ cations there is a propensity for AA water. This creates rings in which Na+ is a double donor and one of the water molecules is a double acceptor. This scenario resembles protonated water72 more than bulk water. The AA interaction consists of two weak HBs (due perhaps to ring strain), which together exceed the energy of a simple HB by ΔEAA(p). Interestingly, ΔEAA(p) is also non-monotonic in p, reaching its maximum for n = 4 (Table 1). We ascribe this 1670


Article


the outer shells), more than any other low-energy isomer; hence, they dominate at room temperature and above. These examples suggest that the concept of free vs restricted rotation of DWMs may be useful in assessing the temperaturedependent free energy ordering. 4.3. Infrared Spectra. We have compared the measured spectra for n = 4 and 5 at low (IRPD) and high temperatures (IRMPD) with computed VPT2 anharmonic spectra. Agreement is good, much better than with the harmonic spectra previously calculated. The more accurate calculations are necessary for the identification of the various isomers contributing to the experimental signal. With the current accuracy levels, positive isomer identification could be achieved for the tetramer, with only tentative results for the pentamer. Recall that for n = 4 and 5 there is a single “lowest energy isomer”, 400 and 410, respectively, that is both the global minimum in electronic energy and free energy, over the whole temperature range. However:

Figure 12. Vibrational normal-mode frequencies (in the harmonic approximation) for the various isomers of the hexa-aqua sodium cation cluster, as a function of the mode number. Inset shows an enlargement of the low-frequency tail. Of the 51 normal modes, 18 are intramolecular water modes: 12 water OH stretches (same as the number of H atoms), above 3000 cm−1, and 6 HOH bending modes, all around 1650 cm−1.

• The low-T (100 K) spectrum is not dominated by the lowest energy isomer. For n = 4 the dominant isomer in the spectrum is 310 (2.4 kcal/mol higher electronic energy) whereas for n = 5 there seems to be a mixture of isomers even at low T, composed of 311, 320-2, possibly also 410, and even 230, which is 12 kcal/mol higher in electronic energy. • Except for the lowest energy isomer, none of the isomers identified spectroscopically is predicted from the Boltzmann distribution of isomers, and vice versa. Thus, the experimental sample is remote from equilibrium. • The CN of the spectroscopically observed n = 4 and 5 isomers tends to average around 3 rather than 4. In order to understand how this may come about, we portray in Figure 13 all the isomers suggested from the IR analysis, with arrows indicating single water molecule loss reactions connecting them, and their electronic energy differences according to the stability rules of Table 1. In this way we obtain a “map” of sequential water-loss pathways connecting the various isomers for n = 3−5. Assuming that small clusters are generated from larger ones along (mostly unimolecular) water-loss pathways, one may understand why 310 dominates the low-T tetramer spectrum, rather than the lowest energy isomer, 400. Most of the pathways from n = 5 to 4 lead to 310,

An important contribution to the entropy comes from internal rotations that, when slightly hindered, become soft vibrational modes. Consider the “dangling” water molecules (DWMs) that, aside to a single acceptor bond (to the Na+ or another water molecule), are not further hydrogen-bonded. The rotation of DWMs in the first solvation shell is hindered substantially by the repulsion between the inward-pointing LPs (section 3.2), and more so for larger CNs. However, the rotation of DWMs in the outer (second or third) solvation shells is less hindered and will contribute substantially to the entropy. In comparison, isomers with a high content of AA bonds are rather rigid. Their population will decline swiftly with increasing T, and isomers with DWMs will take over. Consider again the hexamer isomeric distribution in Figure 7. The lowest electronic energy isomer, 420-1, has two AA bonds and no water molecules that are free to rotate. It dominates at low T, but declines rapidly with increasing T. The first isomer on the electronic energy scale that has a second-shell dangling HB is 411, which is only 1.6 kcal/mol above the global minimum isomer. It thus rapidly increases in population and dominates the whole high-T regime. Between 420-1 and 411 in Figure 6, we find 510 and 600. 510 is slightly lower in energy, and its first-shell water molecules are less restricted in their rotational motion than 600, so it populates slightly with increasing T, whereas 600 not at all. At the high-T end we witness 420-5 rise in population over 420-4, although both have similar energies and one DWM in the second shell. However, 420-5 has two DWMs in its first shell, whereas 420-4 has only one. Continuing to the heptamer, we note in Figure 8 two accidentally degenerate isomers, 520-1 and 421-1. The first has an inner-shell DWM, whereas the second has a third-shell DWM. Due to its free internal rotation, we witness in Figure 9 that while the two isomers are equally probable as T→ 0, the population of the 421-1 isomer increases with increasing T, while that of 520-1 decreases. At higher T, the population of isomers with a single DWM decreases, and isomers with multiple DWMs dominate. In particular, isomers 412 and 421-3 have four DWMs each (two in the inner shell and two in

Figure 13. Putative water-loss pathways connecting the isomers for n = 5, 4, and 3 that have been identified from the IR spectra with the aid of the VPT2 anharmonic frequencies. Blue/red-colored isomers are identified only in the low/high-T spectra, whereas green isomers may appear in both. Number on each arrow indicates the electronic energy difference (in kcal/mol) between the two isomers connected by it, as evaluated from the various HB-type energies of Table 1. Red arrows indicate high-energy pathways, where the pathway from 320-2 to 220-2 is too high in energy to merit indication. 1671


Article


This argues that more work is required before full understanding of the mechanism responsible for the non-equilibrium isomer distribution is achieved. Finally, we note that much work has been done on assigning spectral shifts to various types of HBs in pure water (bulk or clusters).74,75 Generally speaking, the stronger the HB, OH···O, the weaker the covalent OH bond, and thus the more redshifted is its stretching frequency. Less work along these lines is available for aqueous sodium ion clusters. In Figure 14 we show a correlation for the anharmonic VPT2 frequencies obtained here for n = 4 and 5 (data in the SI) with the “stability rules” HB strengths from Table 1. The strongest and most red-shifted HBs are from a water molecule that donates a single HB (its other, “dangling” OH bond is not engaged in a HB). The “simple” and AAD HBs belong to this class. Weaker and less red-shifted HBs are from a water molecule that donates two HBs (our aDD and sDD types). The weakest and least redshifted are the OH stretches from two water molecules in the first shell donating to the same second-shell water (AA-type). Additionally, for the same HB-type, smaller CN leads to somewhat larger HB strengths and red-shifts.

and they involve lower vaporization energies than the other pathways. While a water vaporization model may have the potential for explaining the non-equilibrium isomer distribution, it is not clear whether it is applicable to the experiments of Miller and Lisy.54 In these experiments, neutral water−argon clusters are prepared at very low temperatures with large excess of argon. Na+ prepared by thermionic emission collides with the neutral cluster with excess energy of about 230 kcal/mol. Miller and Lisy54 have assumed that the daughter cluster cools down by evaporating the most labile species, which is Ar. This model has been recently turned into a quantitative DFT-based MD protocol.56 Due to the high computational cost of including so many Ar atoms, the evaporative cooling steps were mimicked by simply scaling down the velocities. Thus, we can discuss the cluster-cooling mechanism only qualitatively at this stage. “Evaporative cooling” means that the excess energy is used solely to overcome the binding energies of the emitted particles (i.e., their relative kinetic energies are negligible). With a water−Ar (ZPE-corrected) binding energy of 1−2 kcal/mol,54 this may require, say, 150 Ar atoms. But then they cannot all bind to a core of only 3−4 water molecules, and the Ar−Ar interaction (not ZPE-corrected) is only 100 cm−1.73 Then, to dissipate the excess energy more than 500 Ar atoms may be needed. Experimentally, it was not demonstrated that so many Ar atoms are involved. Theoretically, the mechanism may be more complicated than anticipated:

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00038. Cartesian coordinates and thermodynamic functions for all of the Na+(H2O)n=1−8 isomers studied herein; harmonic and anharmonic frequencies and intensities for n = 4 and 5; free energy differences for the n = 3 and 4 isomers as a function of temperature (PDF)

(a) Ar atoms that fly out immediately on impact might carry a large excess of kinetic energy, reducing the number of subsequent vaporization events needed for dissipating all of the excess energy. (b) Water molecules may also participate in the cooling process.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the Israel Science Foundation (grant no. 766/12). The Fritz Haber Center is supported by the Minerva Gesellschaft für die Forschung, München, FRG. We thank James M. Lisy for the experimental spectra.

■

REFERENCES

(1) Pohl, H. R.; Wheeler, J. S.; Murray, H. E. Met. Ions Life Sci. 2013, 13, 29−47. (2) Dudev, T.; Lim, C. J. Am. Chem. Soc. 2010, 132, 2321−2332. (3) Payandeh, J.; Scheuer, T.; Zheng, N.; Catterall, W. A. Nature 2011, 475, 353−359. (4) Caminiti, R.; Licheri, G.; Paschina, G.; Piccaluga, G.; Pinna, G. J. Chem. Phys. 1980, 72, 4522−4528. (5) Chizhik, V. I.; Mikhailov, V. I.; Su, P. C. Theor. Exp. Chem. 1987, 22, 480−483. (6) Kameda, Y.; Sugawara, K.; Usuki, T.; Uemura, O. Bull. Chem. Soc. Jpn. 1998, 71, 2769−2776. (7) Mähler, J.; Persson, I. Inorg. Chem. 2012, 51, 425−438. (8) Pálinkás, G.; Riede, W. O.; Heinzinger, K. Z. Naturforsch., A: Phys. Sci. 1977, A32, 1137−1145. (9) Mezei, M.; Beveridge, D. L. J. Chem. Phys. 1981, 74, 6902−6910. (10) Impey, R. W.; Madden, P. A.; McDonald, I. R. J. Phys. Chem. 1983, 87, 5071−5083.

Figure 14. Correlation between OH frequency and OH···O HB strength (colored circles). The frequencies are taken from the VPT2 anharmonic calculations for the n = 4 and 5 isomers reported in the SI, and averaged between symmetric and asymmetric modes when relevant. The energies are obtained from Table 1 as follows: AAD (magenta) ΔEAAD for isomer 311; “simple” (red) ΔEw(2) for isomers 220-1 and 230, ΔEw(3) for isomers 320-1 and 230; aDD (orange) ΔEaDD(3) for isomer 320-2; sDD (green) ΔEsDD(2) for isomers 220-2 and 230; AA (blue) [ΔEw(3) + ΔEAA(3)]/2 for isomers 310 and 320-1, [ΔEw(4) + ΔEAA(4)]/2 for isomer 410. The insets depict schematically four of the HB-types. Red circles are O atoms, OH is black line, whereas the absorbing OH is colored like the corresponding data point. 1672


Article

Journal of Chemical Theory and Computation (11) Chandrasekhar, J.; Spellmeyer, D. C.; Jorgensen, W. L. J. Am. Chem. Soc. 1984, 106, 903−910. (12) Heinzinger, K. Pure Appl. Chem. 1985, 57, 1031−1042. (13) Bounds, D. Mol. Phys. 1985, 54, 1335−1355. (14) Rowley, C. N.; Roux, B. J. Chem. Theory Comput. 2012, 8, 3526−3535. (15) Lev, B.; Roux, B.; Noskov, S. Y. J. Chem. Theory Comput. 2013, 9, 4165−4175. (16) Skipper, N. T.; Neilson, G. W. J. Phys.: Condens. Matter 1989, 1, 4141−4154. (17) Mancinelli, R.; Botti, A.; Bruni, F.; Ricci, M. A.; Soper, A. K. J. Phys. Chem. B 2007, 111, 13570−13577. (18) Megyes, T.; Bálint, S.; Grósz, T.; Radnai, T.; Bakó, I. J. Chem. Phys. 2008, 128, 044501. (19) Kulik, H. J.; Marzari, N.; Correa, A. A.; Prendergast, D.; Schwegler, E.; Galli, G. J. Phys. Chem. B 2010, 114, 9594−9601. (20) Bankura, A.; Carnevale, V.; Klein, M. L. J. Chem. Phys. 2013, 138, 014501. (21) Bankura, A.; Carnevale, V.; Klein, M. L. Mol. Phys. 2014, 112, 1448−1456. (22) Patwari, G. N.; Lisy, J. M. J. Chem. Phys. 2003, 118, 8555−8558. (23) Kistenmacher, H.; Popkie, H.; Clementi, E. J. Chem. Phys. 1974, 61, 799−815. (24) Perez, P.; Lee, W. K.; Prohofsky, E. W. J. Chem. Phys. 1983, 79, 388−392. (25) Lybrand, T. P.; Kollman, P. A. J. Chem. Phys. 1985, 83, 2923− 2933. (26) Arbman, M.; Siegbahn, H.; Pettersson, L.; Siegbahn, P. Mol. Phys. 1985, 54, 1149−1160. (27) Probst, M. M. Chem. Phys. Lett. 1987, 137, 229−232. (28) Cieplak, P.; Lybrand, T. P.; Kollman, P. A. J. Chem. Phys. 1987, 86, 6393−6403. (29) Dang, L. X.; Rice, J. E.; Caldwell, J.; Kollman, P. A. J. Am. Chem. Soc. 1991, 113, 2481−2486. (30) Perera, L.; Berkowitz, M. L. J. Chem. Phys. 1991, 95, 1954− 1963. (31) Lin, S.; Jordan, P. C. J. Chem. Phys. 1988, 89, 7492−7501. (32) Glendening, E. D.; Feller, D. J. Phys. Chem. 1995, 99, 3060− 3067. (33) Merrill, G. N.; Webb, S. P.; Bivin, D. B. J. Phys. Chem. A 2003, 107, 386−396. (34) Mhin, B. J.; Kim, J.; Kim, K. S. Chem. Phys. Lett. 1993, 216, 305−308. (35) Kim, J.; Lee, S.; Cho, S. J.; Mhin, B. J.; Kim, K. S. J. Chem. Phys. 1995, 102, 839−849. (36) Bauschlicher, C. W., Jr.; Langhoff, S. R.; Partridge, H.; Rice, J. E.; Komornicki, A. J. Chem. Phys. 1991, 95, 5142−5148. (37) Hashimoto, K.; Morokuma, K. J. Am. Chem. Soc. 1994, 116, 11436−11443. (38) Lee, H. M.; Tarakeshwar, P.; Park, J.; Kołaski, M. R.; Yoon, Y. J.; Yi, H.-B.; Kim, W. Y.; Kim, K. S. J. Phys. Chem. A 2004, 108, 2949− 2958. (39) Rao, J. S.; Dinadayalane, T. C.; Leszczynski, J.; Sastry, G. N. J. Phys. Chem. A 2008, 112, 12944−12953. (40) Biring, S. K.; Sharma, R.; Misra, R.; Chaudhury, P. J. Cluster Sci. 2013, 24, 715−737. (41) Neela, Y. I.; Mahadevi, A. S.; Sastry, G. N. Struct. Chem. 2013, 24, 67−79. (42) Calvo, F.; Doye, J. P. K.; Wales, D. J. Chem. Phys. Lett. 2002, 366, 176−183. (43) Temelso, B.; Archer, K. A.; Shields, G. C. J. Phys. Chem. A 2011, 115, 12034−12046. (44) Lv, Z.-L.; Xu, K.; Cheng, Y.; Chen, X.-R.; Cai, L.-C. J. Chem. Phys. 2014, 141, 054309. (45) Christie, R. A.; Jordan, K. D. J. Phys. Chem. B 2002, 106, 8376− 8381. (46) Kuo, J.-L.; Klein, M. L. J. Chem. Phys. 2005, 122, 024516. (47) Luo, Y.; Maeda, S.; Ohno, K. J. Phys. Chem. A 2007, 111, 10732−10737.

(48) Nguyen, Q. C.; Ong, Y.-S.; Kuo, J.-L. J. Chem. Theory Comput. 2009, 5, 2629−2639. (49) Temelso, B.; Köddermann, T.; Kirschner, K. N.; Klein, K.; Shields, G. C. Comput. Theor. Chem. 2013, 1021, 240−248. (50) Akase, D.; Teramae, H.; Aida, M. Chem. Phys. Lett. 2015, 618, 51−56. (51) Fifen, J. J.; Nsangou, M.; Dhaouadi, Z.; Motapon, O.; Jaidane, N.-E. J. Chem. Theory Comput. 2013, 9, 1173−1181. (52) Fifen, J. J.; Nsangou, M.; Dhaouadi, Z.; Motapon, O.; Jaidane, N.-E. J. Chem. Phys. 2013, 138, 184301. (53) Miller, D. J.; Lisy, J. M. J. Am. Chem. Soc. 2008, 130, 15393− 15404. (54) Miller, D. J.; Lisy, J. M. J. Am. Chem. Soc. 2008, 130, 15381− 15392. (55) Ke, H.; van der Linde, C.; Lisy, J. M. J. Phys. Chem. A 2015, 119, 2037−2051. (56) Brites, V.; Cimas, A.; Spezia, R.; Sieffert, N.; Lisy, J. M.; Gaigeot, M.-P. J. Chem. Theory Comput. 2015, 11, 871−883. (57) Kulig, W.; Agmon, N. J. Phys. Chem. B 2014, 118, 278−286. (58) Kulig, W.; Agmon, N. Phys. Chem. Chem. Phys. 2014, 16, 4933− 4941. (59) Bloino, J.; Barone, V. J. Chem. Phys. 2012, 136, 124108. (60) van Duijneveldt, F. B.; van Duijneveldt-van de Rijdt, J. G. C. M.; van Lenthe, J. H. Chem. Rev. 1994, 94, 1873−1885. (61) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553−566. (62) Frisch, M. J.; et al. Gaussian 09, Revision D.1; Gaussian Inc.: Wallingford, CT, 2009. (63) Džidić, I.; Kebarle, P. J. Phys. Chem. 1970, 74, 1466−1474. (64) Gillespie, R. J. Coord. Chem. Rev. 2008, 252, 1315−1327. (65) Tielrooij, K. J.; Garcia-Araez, N.; Bonn, M.; Bakker, H. J. Science 2010, 328, 1006−1009. (66) Ludwig, R. Angew. Chem., Int. Ed. 2001, 40, 1808−1827. (67) Laing, M. J. Chem. Educ. 1987, 64, 124−128. (68) Agmon, N. Acc. Chem. Res. 2012, 45, 63−73. (69) Kamarchik, E.; Wang, Y.; Bowman, J. M. J. Chem. Phys. 2011, 134, 114311−114319. (70) Ramaniah, L. M.; Bernasconi, M.; Parrinello, M. J. Chem. Phys. 1998, 109, 6839−6843. (71) Xantheas, S. S. Chem. Phys. 2000, 258, 225−231. (72) Hassanali, A.; Giberti, F.; Cuny, J.; Kühne, T. D.; Parrinello, M. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 13723−13728. (73) Colbourn, E. A.; Douglas, A. E. J. Chem. Phys. 1976, 65, 1741− 1745. (74) Lenz, A.; Ojamäe, L. J. Phys. Chem. A 2006, 110, 13388−13393. (75) Auer, B.; Kumar, R.; Schmidt, J. R.; Skinner, J. L. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 14215−14220.

1673


Structure and Spectroscopy of Hydrated Sodium Ions at Different

Recommend Documents