Global Minima of Protonated Water Clusters - American Chemical

Oct 5, 2012 - (MP2)) agree that for protonated clusters, the cage-like morphology ... and prismatic “edge-sharing pentagonal prism” (EPP) cluster,...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/JPCA

Global Minima of Protonated Water Clusters (H2O)20H+ Revisited P. Parkkinen,* S. Riikonen, and L. Halonen Laboratory of Physical Chemistry, Department of Chemistry, University of Helsinki, P.O. Box 55, FI-00014 Helsinki, Finland

ABSTRACT: A graph-theoretical analysis is performed on the (H2O)20 “edge-sharing pentagonal prism” cluster to find proton configurations that yield the lowest cluster total energies. Using the low-energy structures, we create models for protonated (H2O)20H+ clusters that compete energetically with “cage-like” structures proposed earlier in the literature. We perform benchmarking between different theoretical methods and observe significant stabilization of compact versus polyhedral clusters due to long-range electron correlation effects, which make the comparison between different cluster morphologies difficult using density functional theory only. All methodologies that we used (up to the second level of Møller−Plesset perturbation theory (MP2)) agree that for protonated clusters, the cage-like morphology proposed in [Chem. Phys. Lett. 2000, 324, 279−288] is the most stable one. We study in detail several (H2O)20H+ cluster structures and suggest that the energetics in small protonated water clusters is dominated by the competition of open polyhedral structures, favored by the Eigen H9O+4 species, against van der Waals interaction and “ice rules”, which both favor compact structural motifs such as the prismatic particle. We demonstrate this tendency using ab initio calculations for prismatic, dodecahedral, and cage-like clusters, while the “magic number” cluster (H2O)nH+, n = 21, is observed to minimize all competing energy contributions simultaneously. We emphasize the utility of the cage-like cluster as a model template for ice-related atmospheric reactions and benchmark the GPAW density functional theory code for making such calculations.

■

INTRODUCTION

have stimulated a considerable amount of experimental and theoretical research, while one of the most fundamental physical properties of these clusters is the global minimumenergy structure, which is an important prerequisite in understanding other characteristics. According to semiempirical and ab initio calculations, the lowest-energy structures of neutral water clusters (H2O)n with n ≥ 17 tend to be compact, with a maximum number of fourfold-coordinated water molecules and a minimal number of dangling hydrogen bonds.23 Near n = 20, many of the constituent monomers of a cluster for “prismatic” and “cubic” families obey “ice rules”, that is, molecules accept and donate two hydrogen bonds.7,8 The global minimum of neutral water clusters at n = 20 resulting from second-order Møller−Plesset perturbation theory (MP2) ab initio calculations is a compact and prismatic “edge-sharing pentagonal prism” (EPP) cluster,7,8 while the water dodecahedron, which has been a popular model for charged clusters, lies ∼18 kcal/mol (1 kcal/mol = 4.184 kJ/ mol) higher in energy than EPP.8 This has led the authors of ref

Understanding the hydrogen-bonded structures and chemical properties of neutral and protonated water clusters is a timely subject in solvation chemistry, biology, and atmospheric sciences. Solvation of ions in small clusters is expected to resemble the solvation in bulk-phase solutions, while finite clusters are easier to study both experimentally1−6 and computationally.7−15 In the case of biological systems, small protonated water clusters play an important role in intraprotein proton transport.16 In atmospheric chemistry, protonated water clusters are the most abundant ionic species in the Earth’s atmosphere below altitudes of ∼85 km.17 The role of the species in noctilucent cloud formation has been studied in laboratory experiments, and significant lifetimes of clusters containing more than 20 molecules have been observed.18 Many atmospherically appropriate chemical reactions occur on droplets and aerosol surfaces.17,19,20 Water clusters, surfaces, and even a single water monomer are known to enhance some reactions by stabilizing the transition state.20−22 Furthermore, the water clusters, which are sufficiently small to be treated computationally at a high level, can be utilized in benchmarking semiempirical potentials, which can then be used to simulate larger systems. Not surprisingly, the properties of water clusters © 2012 American Chemical Society

Received: August 1, 2012 Revised: October 5, 2012 Published: October 5, 2012 10826

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835

The Journal of Physical Chemistry A

Article

cage-like structure10,14,15). The global minimum, apart from being of fundamental interest, may be used to test the efficiency of algorithms that sample the potential energy surface. For water systems, this task is difficult. In addition to all possible molecular sites and geometries, that is, polyhedral, cubic, prismatic, and so forth “families”, there is the issue of proton disorder, which makes the energy landscape rugged.7 Proton disorder can also result in energetic overlap between different families of water clusters, which should all be efficiently sampled. In order to make an exhaustive study of the most stable family of water clusters at n = 20 (prismatic family), we have performed a graph-theoretical analysis of the EPP cluster. Earlier, this analysis was performed for the water dodecahedron (i.e., a polyhedral cluster) and for a small cubic water cluster.30 To our knowledge, the present work is the first exhaustive study of proton disorder in a prismatic morphology. A graph-theoretical analysis makes it possible to compare energetics of neutral water clusters that exhibit different morphologies (say, dodecahedral versus prismatic), while keeping each morphology in its proton-configurational energy minimum. This is important as we observe a strong energy dependence on the proton configuration. One can draw artificial conclusions on the relative stability of different morphologies if an arbitrary proton configuration is chosen, a problem that is present in the literature. We further construct protonated water clusters based on optimal proton-configurational neutral clusters and discuss how the cage-like structures found by Hodges and Wales10 combine in an optimal way a tetrahedrally coordinated hexagonal ice-like structure with the quasi-planar H9O+4 Eigen species. This species is seen to be responsable for favoring polyhedral structures and in driving the morphologies of charged clusters into something different from neutral clusters. Finally, we emphasize the possibility of using EPP-type clusters as templates for studying ice-related atmospheric reactions. These clusters are sufficiently large to exhibit at the same time both molecules that obey strictly the ice rules and molecules that have either oxygen or hydrogen “dangling bonds”. They are also small enough to be handled up to the MP2 level of calculation by modern quantum chemistry codes. We have performed benchmarking and comparative calculations of neutral and protonated water clusters using the GPAW DFT31,32 and the TurboMole quantum chemistry program packages,33 with different basis sets. Van der Waals (vdW) forces are observed to play an important role in stabilizing compact water cluster morphologies.

8 to suggest in the favor of the prismatic cluster that the neutral water dodecahedron is unlikely to be abundant under atmospheric conditions.8 In experiments, which are typically performed on protonated (H2O)nH+ clusters, a different picture has emerged; the first mass distribution spectra, published as early as 1973,24 featured a sharp peak at n = 21, and later on, an open structure, a water dodecahedron enclosing a single molecule, was proposed.25,26 The water dodecahedron or a similar (broken) polyhedral structure with three-fold-coordinated water molecules enclosing a single water monomer or a hydronium ion (i.e., a clathrate structure, analogous to natural gas hydrates27) has since been discussed in innumerous experimental2−5 (see also ref 1 and references therein) and theoretical9−15 works. A dodecahedral clathrate structure with three-fold-coordinated water molecules was appealing because it can account for key experimental features of charged clusters; infrared-spectroscopy (IR) studies (refs 2−5 and ref 1 and references therein) show fingerprints of three-fold water molecules at n = 21. In more detail, the IR spectra as a function of the cluster size show a “collapse” of spectral features into a single peak at n = 21, which is associated with the OH stretches of three-foldcoordinated water molecules.2−4 The number of dangling hydrogen bonds has also been tested using a trimethylamine attachment to the cluster, which has shown a clear peak with 21 water molecules and 10 dangling hydrogen bonds. The latter number fits into the dodecahedron picture.28 Several theoretical calculations on global minima structures of protonated water clusters (H2O)nH+ at n = 20 and 21 now exist,2,9−15,29 the used methodology including density functional theory (DFT)9 and MP22,11,14,15 ab initio calculations, KJ (Kozack−Jordan) and ASP (asymmetric site) potentials,10 the OSS2 potential,13 and the EVB (empirical valence bond) method.12 Most of these studies,9,10,12,14,15 where the global minimum structure at n = 20−21 is found with a structural search9,10,12 or is taken from the work of Hodges and Wales,10 share common features, namely, they contain one H3O+ as part of the surface and a polyhedral cage-like structure, which binds to a central water monomer. For this reason, these clusters can be at the same time characterized as “polyhedral” (the outer cluster surface) or “compact” and “bulk-like” as the central water monomer binds to the cluster surface, fullfilling ice rules, that is, donating and accepting two hydrogen bonds. These clusters typically have nine dangling hydrogen bonds projecting into the vacuum, and at n = 21, a perfect distorted water dodecahedron with a central water monomer emerges.10 In the following, we shall refer to these structures as “cage-like”. Protonated clusters at finite temperatures have been explored using the OSS213 potential combined with basin-hopping and parallel-tempering algorithms. At n = 20 and 21, distorted polyhedral cages typically emerge at low temperatures, higher temperatures favoring more open structures, “trees” and “flowers”.13 Picosecond-scale ab initio molecular dynamics at 150 and 300 K have been performed by Bankura et al. at n = 21, and an increase in low-coordinated water molecules was observed with increasing temperature.15 As emphasized in ref 10, the global minima of protonated cage-like water clusters appear to be different from those of neutral prismatic clusters. This is an interesting point, which we shall explore in detail in this work, taking as a model system the cluster containing 20 water monomers. In the literature, there are inconsistencies concerning the global minimum-energy structure of (H2O)nH+ at n = 20 (EPP,11 dodecahedron,13 or

■

RESULTS Edge-Sharing Pentagonal Prism. Some EPP water clusters7,8 having different hydrogen bond topologies are depicted in Figure 1. Taking a closer look at Figure 1, the name of the cluster becomes obvious. It consists of three prisms, each one having facets of five water molecules and edges of four water molecules. Two edges in each prism are shared, while the remaining edges face the vacuum, exhibiting a varying number of different dangling hydrogen bonds. The water molecules at the edges can be classified as follows: molecules having a dangling hydrogen atom are double acceptor, single donor (AAD) species, while molecules exhibiting a dangling oxygen atom are denoted as single acceptor, double donor (ADD) species. Molecules in the “interior” of the cluster (eight molecules in total) obey the ice 10827

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835

The Journal of Physical Chemistry A

Article

related to the proton order versus proton disorder of a cluster, a zero dipole moment indicating proton-disordered clusters. Taking a close look at Figure 1d, it is obvious that large total dipole moments and proton ordering are related to the concentration of (AAD) or (ADD) species on either side of the cluster. The dipole moment distribution of the final sample of 1282 water molecules is presented in Figure 2a. Interestingly, the

Figure 1. Some hydrogen bond topologies of the EPP cluster (a−d) and the water dodecahedron (e). In panels (b−e), water molecules containing dangling oxygen bonds (ADD species) and dangling hydrogen bonds (AAD species) are marked with blue and green, respectively. Otherwise, oxygen and hydrogen atoms are marked with red and white, respectively. Black lines denote hydrogen bonds. (b) A low-energy hydrogen bond topology resulting from our graphtheoretical search (with a dipole moment of 0.8 D). (c) A low-energy hydrogen bond topology exhibiting a larger dipole moment of 4.1 D. (d) Hydrogen bond topology with an extreme dipole moment of 24.5 D.

Figure 2. (a) Dipole moment distribution of the sample and the Boltzmann distribution weighted dipole moment distributions at 203 (b) and 273 K (c).

sample contains a large number of clusters exhibiting huge dipole moments, the largest peak occurring at ∼9 D. After weighting the distribution with Boltzmann weights using total energies from PBE DFT calculations, only two nominal peaks remain, a main peak at ∼0.0−1.0 D and a smaller peak at ∼4.0−5.0 D. According to Figure 2c, at T = 273 K, approximately 5% of the clusters exhibit a slight proton ordering and a considerable dipole moment of ∼4.5 D. Zero dipole moment (proton-disordered) clusters have increased stability, as is seen in Figure 3a, where the cluster total energy as a function of the total dipole moment is plotted. The plot is scattered, but a trace of quadratic dependence can still be observed, while a clear quadratic dependence has been observed for the water dodecahedron for which the clusters with the most favorable stabilization energies also have the lowest total dipole moments.36,37 The group of cluster topologies with M ≈ 0 D in Figure 3a (indicated with a dashed box) includes 101 configurations of the n = 1282 sample, while for the true n = 22 960 sample, we estimate this number to be about 950 configurations as approximated by summing up the dipole moments of individual water monomers. A few representative cluster geometries from the scatter plot of Figure 3a have been visualized in Figure 1; the low-energy structures near M = 0 D and 5 D can be found in panels (b) and (c), respectively, while the highest-energy structure at M ≈ 25 D is depicted in panel (d). From the sample of n = 1282 cluster geometries, we chose a reduced sample of n = 8 representative geometries that span the values of dipole moments from ∼0 D to the maximum of ∼25 D. The objective of this smaller sample is to provide cluster geometries for benchmarking ab initio and quantum chemistry codes. Figure 3b,c shows dipole moments and relative stabilization energies per molecule of the n = 8 sample set computed with GPAW (with the PBE exchange and correlation functional) and TurboMole using DFT with PBE or a B3LYP exchange and correlation functional employing two different

rules, that is, they both accept and donate two water molecules (AADD). There are six (AAD) and (ADD) species in the EPP cluster, which can be distributed, by choosing an adequate proton configuration, on different edges in an arbitrary manner. Clusters in Figure 1 exhibit different dangling bond distributions; in panels (b) and (c), the two (AAD) species on the rightmost edge of the cluster have different arrangements. Panel (d) features an extreme case where the lowermost edge has a high concentration of (ADD) species (dangling oxygen atoms) and the other edges have a high concentration of (AAD) species. A graph-theoretical analysis on the hydrogen bond topology of the EPP cluster was performed following ref 34. Altogether, 22 960 symmetrically distinct proton configurations were found, from which 1282 cluster geometries were randomly chosen and optimized using the GPAW DFT code. The most optimal hydrogen bond topology found this way, EPP1, is depicted in Figure 1b. It is worth noting that the total number of isomers in the present case is smaller than that for the water dodecahedron (30 026 in ref 30). A smaller number of isomers for EPP results from the larger number of hydrogen bonds in this cluster, leading to more constraints. On the other hand, the number of isomers for EPP is increased due to the lower symmetry (D3h) when compared to the dodecahedron (Ih), the overall net result being fewer isomers for EPP. Apart from the total energy, a key physical quantity in our analysis is the total dipole moment. For water clusters, the total dipole moment can be given as an estimate by summing up the dipole moments (M) of individual water molecules. In water, M = 1.855 D (1 D = 3.33564 × 10−30 C m), while in clusters, such as the prismatic hexamer, M ≈ 2.7 D.35 The dipole moment is 10828

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835

The Journal of Physical Chemistry A

Article

Figure 3. (a) Stabilization energies (eq 1) per water monomer (kcal/mol) as a function of the dipole moment of the EPP cluster in a sample of 1282 structures. Panel (b) presents dipole moments and panel (c) stabilization energies per water monomer (kcal/mol) of eight sample clusters calculated using different schemes. In panel (c), the stabilization energies per molecule have been shifted by E1 (=stabilization energy per molecule of the sample number 1). In panels (b) and (c), the GPAW calculations with PBE are indicated by a blue line (E1 = −10.70 kcal/mol), while the rest of the calculations have been performed using TurboMole as follows: PBE and the def2-TZVPD basis set (a green line, E1 = −10.66 kcal/mol), B3LYP and the def2-TZVPD basis set (a red line, E1 = −9.33 kcal/mol), and B3LYP and the def2-SVPD basis set (a light blue line, E1 = −11.91 kcal/mol).

Table 1. Relative Stabilization Energies of Some Neutral (eq 1) and Protonated (eq 2) (H2O)20 Clusters (see Figures 1b and 4) as Calculated with GPAW (PBE Functional) and TurboMole (B3LYP Functional or MP2 with Either the def2-TZVPD or the def2-SVPD Basis Setsa method (basis set) structure

GPAW/PBE

B3LYP (SVPD)

B3LYP (TZVPD)

MP2 (SVPD)

MP2 (TZVPD)

EPP1 EPP2 EPP1(+) EPP1(+)† EPP2(+) EPP3(+) WH(+) D D(+)

76.39 76.64 9.27 9.19 5.91 7.02 0.00 (−287.49) 77.13 6.30

77.78 77.71 11.55 9.73 8.76 9.45 0.00 (−287.57) 82.73 13.46

76.48 76.08 (*) 8.52 7.68 7.82 0.00 (−263.04) 77.54 9.73

77.44 78.94 13.37 13.80 10.63 12.55 0.00 (−328.32) 91.91 24.67

77.14 77.77 (*) 12.23 9.77 10.32 0.00 (−295.48) 85.51 18.18

a

Absolute values of stabilization energies are indicated in parentheses. All energies are in kcal/mol. Structure WH(+) has been taken from refs 10 and 48, while D is from ref 30. Geometries marked with (*) are unstable. For the meaning of EPP1(+)†, see the text.

In previous studies, there has been some controversy regarding the energy ordering of protonated clusters. According to Khan’s studies,11 the zero-temperature energy ordering of clusters is maintained once these clusters are protonated as the (H2O)20H+ dodecahedron lies 22.0 kcal/mol higher in energy than protonated EPP. However, the MP2 calculations of Kus et al.14 showed that a protonated cage geometry proposed in ref 10 was ∼4 kcal/mol more favorable than the protonated EPP cluster. Kus et al. also suggested that the energy difference between the protonated EPP and dodecahedron structures is small (∼4 kcal/mol) and that the results were sensitive with respect to the basis set size and the method used (HF versus MP2). Studies employing semiempirical potentials identified the (H2O)20H+ global minimum as either a dodecahedron13 or a cage-like geometry.10 In ref 10, the authors emphasize the absence of prismatic structures similar to those found by Khan.11 The low-energy proton configuration EPP1 (Figure 1b) presented in the previous section is now employed. Analysis of the proton disorder is mandatory to make a consistent comparison against other cluster morphologies. According to Figure 3a, taking a high-energy proton configuration within low dipole moment clusters may create an energy “jump” of more than 8 kcal/mol (=20 × 0.4 kcal/mol) for a cluster containing

basis sets. From Figure 3b,c, we observe that if one is interested in the relative energies and energy ordering of clusters, all schemes give consistent results. Furthermore, the curve of Figure 3b can be recovered by a vector sum over monomer dipole moments and by scaling the monomer dipole moments to 3.93 D. This implies that the first screening of low-energy structures may be done using the classical vector sum of dipoles. Finally, according to Figure 3, there is no significant error associated with solving the Poisson equation with zeroboundary conditions (a technique employed by GPAW) for clusters exhibiting large dipole moments. Global Minima: EPP, Dodecahedral, and Cage-like. Among the structures of neutral (H2O)20 clusters, the superior stability of the EPP cluster at 0 K with respect to other morphologies such as cubic, dodecahedral, and so forth is wellestablished in the literature,7,8 and compact structures in general, such as prismatic and cubic, are more stable than open structures. Fanourgakis et al.8 found the energy ordering with MP2 calculations for prismatic, cubic, and dodecahedral geometries to be 10.90, 10.63, and 10.01 kcal/mol per molecule, respectively. The (H2O)20 dodecahedron was found to lie 17.8 kcal/mol higher in energy than EPP. 10829

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835

The Journal of Physical Chemistry A

Article

Figure 4. (a−c) Creating optimal (H2O)20H+ geometries from the minimum-energy structure of the neutral (H2O)20 cluster. (a) The minimumenergy proton-ordering configuration of the EPP cluster (see also Figure 1a) with the extra proton inserted into a dangling oxygen site. (b) The structure of panel (a) modified; the molecule marked with a dashed arrow has been rotated 180°, and the extra proton has moved to a new site above the rotated molecule. In panels (a) and (b), the solid arrows mark the start and end points of spontaneous proton migration during geometry optimization. (c) The structure of panel (b) modified by proton migration and by a slight rearrangement of the proton ordering. (d) The cage-like cluster structure proposed by Hodges and Wales;10,48 ADD (blue), AAD (green), AD (black), and H3O+ (yellow) species have been marked with different colors. In panel (d), a simplified geometry is also shown where only the tetrahedrically coordinated species and the Eigen (H9O+4 ) complex are visible. (e) A protonated water dodecahedron constructed by adding a proton to the optimal neutral cluster from ref 30.

via the Grothuss mechanism,39,40 the Zundel H2O+5 form acting as an intermediate state. Nuclear quantum effects (NQEs) blur the distinction between Eigen and Zundel forms as the hydrogen atom becomes delocalized between neighboring oxygens (for more details, see refs 41 and 42 and references therein). The EPP(+) clusters in Figure 4 have been constructed on the basis of the neutral EPP1 cluster (Figure 1b) by inserting a proton into an available AAD species, which transforms it to a DDD species, and then migrating this DDD species by the Grothuss mechanism (either spontaneous migration during geometry relaxation or by rearranging the protons “manually”). For EPP2(+) and EPP3(+), the ice rules have been locally broken by rotating a single monomer, this way creating a site for a DDD species. The geometries of panels (a−c) share a common feature; the overall compact and prismatic structure of the EPP cluster is maintained, while the H3O+ species breaks the ice rules only locally. This way, the energetically favorable ice rules are maintained simultaneously with the triple donor configuration requirement for the H3O+ species. In addition to the triple donor configuration for the H3O+ hydronium species and in thinking in terms of this species only, the energetics of the surrounding Eigen H9O+4 complex as a whole plays an important role in the cluster energetics. The first DFT study of the hydrated proton43 has documented delocalized electronic states (i.e., not localized simply on the H3O+ species), while recently, considerable charge transfer between the H9O+4 complex and its surroundings has been deduced from theoretical calculations.44 In the calculations of protonated prismatic clusters, Khan11 has noted modified O−O distances in the H9O+4 structure. In Table 2, we have listed energies of isolated Eigen structures from several cluster

20 molecules, while incorrectly choosing a cluster with a high dipole moment can make the energy more than 16 kcal/mol higher when compared to the lowest-energy proton configuration within the same morphology. In Table 1, we have compared neutral prismatic (EPP1) and dodecahedral (D) clusters, setting both of them in low-energy proton configurations. At the MP2 level, the prismatic particle is more stable than the dodecahedron, but by only ∼8 kcal/mol, slightly less than that in ref 8, where a more pronounced energy separation of 17.8 kcal/mol was given. We emphasize that comparing energy values between different papers in the literature is not straightforward as different proton configurations are used and, as mentioned earlier, arbitrary and considerable energy “jumps” can be created as a function of the proton disorder. We have used the best known minimumenergy proton configurations for both the neutral dodecahedron (McDonald et al. in ref 30) and the neutral prismatic particle as found by us. Hydrated Proton. When constructing protonated EPP clusters based on the low-energy proton configurations, a few facts about the hydrated proton must be kept in mind; the EPP cluster of Figure 1 is an example how water systems prefer fourfold-coordinated species with a maximum number of molecules that simultaneously accept and donate two hydrogen bonds (AADD species obeying the ice rules). On the contrary, the H3O+ species is a pure hydrogen bond donor (triple donor or DDD). It achieves the lowest-energy state when surrounded by hydrogen bond acceptors only, that is, when forming a H9O+4 Eigen complex.38 These different preferences of the water molecule and the H3O+ ion are known to create the ratelimiting step of proton migration in liquid water; a protontransfer event is preceded by the breaking of a hydrogen bond in the target molecule where the extra proton is about to jump 10830

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835

The Journal of Physical Chemistry A

Article

When the particle is protonated, as described above, the quasi-planar H9O+4 structure should be accommodated into it, which is not simple taking into account the sharp corners accommodating the ADD species. It is then not surprising that according to Tables 1 and 2, the energy ordering of protonated prismatic clusters results from the energy ordering of their H9O+4 species. The local stability of the EPP-based protonated structures was further tested by running low-temperature (200 K) ab initio molecular dynamics using GPAW over 0.5 ps. It was noticed that EPP2(+) and EPP3(+) maintained their geometries, while in EPP1(+), the Eigen species broke some of the hydrogen bonds with the neighboring water molecules, resulting in a slightly modified cluster structure EPP1(+)†. Looking at Table 2, we observe that this behavior is yet another manifestation of the tendency of the Eigen complex to prefer quasi-planar structures and to modify the surrounding cluster accordingly. Both the energy and the parameter d obtain lower values for the modified EPP1(+)†. Dodecahedron. The water dodecahedron (see Figures 1e and 4e) does not possess even a single molecule that would obey the favorable ice rules, and all of its constituent monomers are three-fold-coordinated (i.e., AAD and ADD species). On the other hand, the O−H−O angles are near-perfect, as suggested by Table 2. When the dodecahedron is protonated, it offers, according to Table 2, an ideal site with an optimal H9O+4 structure. Hodges and Wales Structure. The cage-like structure found by Hodges and Wales (WH(+))10 is depicted in Figure 4d. It can be described as a particle having an Ice-Ih interior as there are five molecules that obey the ice rules and are in a nearperfect tetrahedral coordination; the O−H−O angles between

geometries. A tendency showing preference for quasi-planar geometries can be seen. Table 2. Features of Protonated Clusters (GPAW/PBE Calculations)a

Panel (a) defines the “flatness” parameter d for local H9O+4 structures, and panel (b) illustrates the definition of the O−H−O hydrogenbonding angle α. In the table, the mean deviation from α = 180° and the number of ice-rule-obeying monomers (N) (i.e., AADD species) are listed. In the last column, results for the energetics (E) and the height of the H3O+ species from the plane defined by three nearest neighbors (d) of each structure are given. 1 Å = 10−10 m. a

Edge-Sharing Pentagonal Prism. As already mentioned, there are several molecules on the EPP particle that obey the ice rules. However, the geometry of the EPP particle is such that optimal hydrogen bonds are difficult to achieve, as suggested by Khan.11 One measure for the optimization of the hydrogen bond is its O−H−O angle (see panel (b) in Table 2), while stronger hydrogen bonds in the water systems have stronger preference for α = 180°.45 The central molecules in the EPP particles, although obeying the ice rules, deviate up to 17° from the optimal angle of α = 180°.

Figure 5. Results of Table 1 visualized, that is, relative stabilization energies have been plotted for different structures using different methods. In the left inset, each combination of method and basis set is indicated by a symbol/color combination. For each method (and basis set), the reference (i.e., zero-level) energy is the energy of the WH(+) structure. In the right inset, the energy contributions from the Hartree−Fock (Escf) and many-body perturbation (Embp) parts for the neutral clusters have been plotted (MP2 calculations with TZVPD) using EPP1 as the reference system. 10831

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835

The Journal of Physical Chemistry A

Article

the central molecules deviate only 0.8−4.6° from α = 180°. IceIh is an infinite, three-dimensional structure. If one tries to create a finite cluster with it, a termination is needed. For protonated clusters, the quasi-planar Eigen H9O+4 complex is a good candidate. Both Figure 4d and Table 2 suggest that the WH(+) structure combines in an ideal way the energetically optimal three-dimensional Ice-Ih with the quasi-planar Eigen complex. Global Energy Minima. Energetics for several cluster geometries are in Table 1, and the results have been visualized in Figure 5. Different methods/basis sets inequivocally agree that the WH(+) structure is, within neutral and protonated clusters, the most stable structure by at least ∼6 kcal/mol. We attribute this to its ability to combine the hexagonal ice structure with a planar Eigen complex, as described above. For this reason, within each method and basis set, the WH(+) structure has been taken as the reference system. Looking at Figure 5, we observe that when it comes to the relative energy ordering between different systems, all methods give consistent results, that is, WH(+) is energetically the lowest geometry, followed by EPP2(+), then EPP3(+), and so forth. However, the energetic spacing between different structures is sensitive to the basis set and method used. The smaller basis set (def2-SVPD) exaggerates energy differences, while MP2 calculations give systematically a broader energy separation than DFT with PBE and B3LYP. The water dodecahedron is a problematic case as it is sensitive to both the method and the basis set used. However, when taking as the reference benchmark the MP2 calculations using the def2TZVPD basis set, MP2/def2-SVPD is good for the correct energy ordering. To investigate the origin of the significant energy difference (∼8 kcal/mol) between different methods for the neutral dodecahedron (D), we have decomposed its SCF-MP2 total energy into the SCF (Hartree−Fock) and many-body perturbation (MP2) part in the inset of Figure 5, where a clear jump in the MP2 correlation energies for the dodecahedron is observed. The situation looks similar to the water hexamer case, where low-coordinated isomers become unstable once the long-range correlation effects, that is, the vdW interactions, are taken into account.46,47 The vdW interaction stabilizes compact geometries, with a downward shift in energy, an effect that can be qualitatively described by summing up −C6/r6 energy terms over atomic pairs.46 There is a tendency to collapse the water dodecahedron due to the vdW forces. According to Figure 5, the energy spacings, when comparing the PBE and B3LYP calculations, are more pronounced for charged systems. This indicates that the self-interaction error is larger in the case of protonated clusters. However, this is not critical, at least for the compact systems. Taking as a benchmark the MP2/TZVPD calculations, the energy spacing between WH(+) and a protonated prismatic cluster (EPP2(+)) becomes 9.77 kcal/mol (reported as 3.8 kcal/mol in ref 14), while the energy spacing between the protonated dodecahedron and prismatic particle is 8.41 kcal/ mol (3.9 kcal/mol in ref 14 and 22.0 kcal/mol in ref 11). We emphasize that we have used a higher-quality basis set. The 631(1)G*(*) basis sets used in refs 11 and 14 lack additional diffuse orbitals, that is, “aug-” in the Dunning and “D” in the TurboMole notation. In the present work, we have also used systematically proton configurations that give the lowest total energies. For a thorough study of relative energies and energy

overlap between different morphologies, a complete graphtheoretical analysis should be performed in all cases.

■

CONCLUSIONS Sampling the potential energy surface (PES) of water clusters is far from trivial. As pointed out by Wales and Hodges,7 simultaneous sampling of both the center of mass (i.e., oxygen location) and angular (i.e., proton configuration) degrees of freedom is a challenging task as the proton disorder makes the energy landscape rugged. Therefore, low-energy minima within similar morphologies are not necessarily connected by lowenergy pathways. In the present work, a fine-tuned global minimum for the edge-sharing pentagonal prism water cluster family, which is known as the global minima family for (H2O)20 neutral clusters, was searched with the aid of graph-theoretical techniques, that is, by exploring all possible proton configurations. On the basis of this new global minimum of neutral clusters, we further constructed global minima for protonated (H2O)20H+ clusters. Earlier results concerning the sensitive energetics of different morphologies (prismatic, dodecahedral, and cage-like) were recalculated and clarified. All methods considered as well as the present and earlier MP2 calculations agreed that the energy spacing between global minima of different morphologies appeared small, ∼4−10 kcal/mol. We also analyzed why the global minima of neutral and protonated water clusters appear different, as emphasized in ref 10, and suggested that cluster total energy and stability are dominated by the local geometry of the Eigen H9O+4 complex (not simply by the geometry of the H3O+ species), which in turn is described by a single parameter describing the planar nature of the complex. The Eigen species prefers flat geometries, which are readily available in polyhedral cage-like clusters but not found in the prismatic structures. The cage-like cluster found by Hodges and Wales10 was seen to combine in an optimal way a hexagonal ice-like interior with a quasi-planar H9O+4 surface complex. This cluster structure could also be a good computational model for proton-rich ice/water surfaces that are important in atmospheric chemistry.49 Different computational methods (DFT/PBE, B3LYP, and MP2) were tested in the context of neutral and protonated clusters. Stabilization of compact prismatic structures due to van der Waals forces similar to water hexamers46 was observed, although this aspect should be studied in more detail in the future. In this context, our results suggest that the (H2O)20 clusters could be a testbed for studying the validity of dispersion-corrected density functionals50 for water systems,47,51 as clusters exhibiting various levels of “sparsity” may be constructed using the various morphologies of (H2O)20. When performing ab initio molecular dynamics with DFT for small water clusters,15 dispersion corrections should be included. Without them, open structures, which van der Waals forces tend to collapse into more compact motifs, are incorrectly described. In conclusion, at the zero-temperature global minima (exluding zero-point vibrational contributions) of neutral clusters, open, polyhedral cage structures become unstable when compared to compact and prismatic structures due to van der Waals interactions and low-coordinated water monomers. However, when an extra proton is added to the cluster, the resulting H9O+4 species prefers a quasi-planar local structure, which makes polyhedral cage structures energetically appealing as the surface of the polyhedral structures is planar. At the magic number of 21 water molecules, when a perfect 10832

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835

The Journal of Physical Chemistry A

Article

(distorted) polyhedral is formed, various energy-terms obtain simultaneously their minima as the structure is compact (van der Waals contribution), planar on the surface (H9O+4 species) and the central molecule and its neighbours are four-foldcoordinated (fullfilling the ice rules).

Table 3. Water Monomer Dipole Moment Calculated with Different Methods and Basis Setsa

■

METHODS Structural relaxations and molecular dynamics (MD) simulations were performed using the Atomic Simulation Environment52 (ASE, version 3.4). The chosen structural relaxation scheme was BFGSLineSearch. A Verlet algorithm with a time step of 0.2 fs was used for the MD time integration. All MD simulations were performed in the microcanonical (constant total energy) ensemble. Two programs were utilized for the electronic structure part. Most of the DFT calculations were performed with the GPAW code (version 0.8),31,32 while some were repeated with the TurboMole program package.33 The grid spacing of the GPAW calculations was 0.2 Å, and there was 5 Å of vacuum around the simulated cluster. For charged systems, a small Fermi smearing of 50 meV (1 eV = 96.49 kJ/ mol) was introduced in the occupation numbers to converge the electronic structure in addition to the standard Pulay mixing scheme. All B3LYP53 and second-order Møller−Plesset perturbation theory (MP2)54 calculations were done with TurboMole. In these computations, two different basis sets were used, namely, the def2-TZVPD (triple-ζ valence with polarization and diffuse functions) and the def2-SVPD (split valence with polarization and diffuse functions) bases. These basis sets are based on contracted Gaussian-type orbitals (CGTOs),55 and they include diffuse functions56 (D), which are vital in describing static polarizabilities of atoms and molecules correctly. This is necessary for simulating water systems. Our “benchmark” basis set for the MP2 calculations is def2-TZVPD. For neutral water clusters in the MP2 calculations, even a double-ζ basis set with diffuse functions (aug-cc-pVDZ) has been observed to give correct relative energies.8 For all DFT calculations, the dscf module57 was used. The same module was also utilized for the Hartree−Fock part of the MP2 calculations, while for calculating the MP2 correlation energies and forces, the ricc2 module58,59 was employed. The ricc2 module takes advantage of the powerful resolution of the identity (RI) approximation for calculating Coulombic integrals. For def2-SVPD, we used def2-SVP as the RI auxiliary basis set as def2-SVPD RI auxiliary basis sets were unavailable. In ref 56, it was observed that leaving out diffuse basis functions from the auxiliary basis did not affect the overall results, at least in the DFT calculations. In the context of water and ice systems and DFT calculations, a hoard of different gradient-corrected density functionals, combining different exchange and correlation parts, have been used (see, for example, ref 40 and references therein). We have chosen the PBE functional60 (which is similar to PW9161) as it gives excellent results for hexagonal ice62 and water systems.39 As a test, we have calculated the dipole moment of a water monomer, using the different schemes described in this section. These results are listed in Table 3. We observe identical trends as in ref 63, that is, with the PBE functional that we obtain almost an identical dipole moment of ∼1.82 D (1.81 D with BP86 and BLYP functionals in ref 63), while MP2 yields a slightly higher dipole moment of ∼1.88 D (1.88 D in ref 63 with the aug-cc-pVDZ basis set). We obtain excellent

method (basis set)

μ (D)

relative error (%)

GPAW/PBE B3LYP (SVPD) B3LYP (TZVPD) MP2 (SVPD) MP2 (TZVPD) exp. (ref 64)

1.818 1.881 1.851 1.902 1.879 1.855

−1.9 1.5 −0.2 2.6 1.3

a

See the text for details. Relative errors with respect to the experimental value are shown.

agreement with the experimental value using the B3LYP functional and the B3LYP (TZVPD) basis set, while the smaller basis set (SVPD) yields a higher dipole moment. Water clusters in large proton configuration samples of 1282 clusters were relaxed using an acceptance force tolerance of 0.05 eV/Å in order to speed up the calculations. In all other cases (i.e., in smaller samples and in all B3LYP and MP2 calculations), a more accurate force tolerance of 0.02 eV/Å was required before declaring a certain geometry as “stable”. Stabilization energies Es of neutral clusters (H2O)n were calculated as Es = E[(H 2O)n ] − nE[H 2O]

(1)

while for charged clusters (H2O)nH+, the energies were computed as Es = E[(H 2O)n H+] − ((n − 1)E[H 2O] + E[H3O+]) (2)

where E[X] is the total energy of the species X. Zero-point vibrational energy corrections were not taken into account because according to ref 14, the effect of these corrections is only ∼1−2 kcal/mol on relative energy ordering for clusters of 20−21 water monomers.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: pauli.parkkinen@helsinki.fi. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank the Center for Scientific Computing (CSC) for the use of its computational resources. This work has been supported by the Academy of Finland through the FiDiPro, the Center of Excellence (2006−2011), and Lastu programs. We acknowledge useful discussions with Andris Gulans, Kari Laasonen, Nino Runeberg, and Benny Gerber.

■

REFERENCES

(1) Nagashima, U.; Shinohara, H.; Nishi, N.; Tanaka, H. Enhanced Stability of Ion−Clathrate Structures for Magic Number Water Clusters. J. Chem. Phys. 1986, 84, 209−214. (2) Shin, J.-W.; Hammer, N. I.; Diken, E. G.; Johnson, M. A.; Walters, R. S.; Jaeger, T. D.; Duncan, M. A.; Christie, R. A.; Jordan, K. D. Infrared Signature of Structures Associated with the H+(H2O)n (n = 6 to 27) Clusters. Science 2004, 304, 1137−1140. (3) Miyazaki, M.; Fujii, A.; Ebata, T.; Mikami, N. Infrared Spectroscopic Evidence for Protonated Water Clusters Forming Nanoscale Cages. Science 2004, 304, 1134−1137. 10833

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835

The Journal of Physical Chemistry A

Article

(4) Zwier, T. S. The Structure of Protonated Water Clusters. Science 2004, 304, 1119−1120. (5) Douberly, G. E.; Ricks, A. M.; Duncan, M. A.Infrared Spectroscopy of Perdeuterated Protonated Water Clusters in the Vicinity of the Clathrate Cage. J. Phys. Chem. A 2009, 113, 8449−8453. (6) Gutberlet, A.; Schwaab, G.; Birer, O.; Masia, M.; Kaczmarek, A.; Forbert, H.; Havenith, M.; Marx, D. Aggregation-Induced Dissociation of HCl(H2O)4 Below 1 K: The Smallest Droplet of Acid. Science 2009, 324, 1545−1548. (7) Wales, D. J.; Hodges, M. P. Global Minima of Water Clusters (H2O)n, n ≤ 21, Described by an Empirical Potential. Chem. Phys. Lett. 1998, 286, 65−72. (8) Fanourgakis, G. S.; Aprà, E.; Xantheas, S. S. High-Level Ab Initio Calculations for the Four Lowlying Families of Minima of (H2O)20. I. Estimates of MP2/CBS Binding Energies and Comparison with Empirical Potentials. J. Chem. Phys. 2004, 121, 2655−2663. (9) Laasonen, K.; Klein, M. L. Structural Study of (H2O)20 and (H2O)21H+ Using Density Functional Methods. J. Phys. Chem. 1994, 98, 10079−10083. (10) Hodges, M. P.; Wales, D. J. Global Minima of Protonated Water Clusters. Chem. Phys. Lett. 2000, 324, 279−288. (11) Arshad; Khan. Ab Initio Studies of (H2O)20H+ and (H2O)21H+ Prismic, Fused Cubic and Dodecahedral Clusters: Can H3O+ Ion Remain in Cage Cavity? Chem. Phys. Lett. 2000, 319, 440−450. (12) James, T.; Wales, D. J. Protonated Water Clusters Described by an Empirical Valence Bond Potential. J. Chem. Phys. 2005, 122, 134306. (13) Kuo, J.-L.; Klein, M. L. Structure of Protonated Water Clusters: Low-Energy Structures and Finite Temperature Behavior. J. Chem. Phys. 2005, 122, 024516. (14) Kus, T.; Lotrich, V. F.; Perera, A.; Bartlett, R. J. An Ab Initio Study of the (H2O)20H+ and (H2O)21H+ Water Clusters. J. Chem. Phys. 2009, 131, 104313. (15) Bankura, A.; Chandra, A. A First Principles Theoretical Study of the Hydration Structure and Dynamics of an Excess Proton in Water Clusters of Varying Size and Temperature. Chem. Phys. 2011, 387, 92− 102. (16) Garczarek, F.; Gerwert, K. Functional Waters in Intraprotein Proton Transfer Monitored by FTIR Difference Spectroscopy. Nature 2006, 439, 109−112. (17) Niedner-Schatteburg, G.; Bondybey, V. E. FT-ICR Studies of Solvation Effects in Ionic Water Cluster Reactions. Chem. Rev. 2000, 100, 4059−4086. (18) Yang, X.; Castleman, A. W., Jr. Laboratory Studies of Large Protonated Water Clusters Under the Conditions of Formation of Noctilucent Clouds in the Summer Mesopause. J. Geophys. Res. 1991, 96, 22573−22578. (19) Kolb, C. E.; Cox, R. A.; Abbatt, J. P. D.; Ammann, M.; Davis, E. J.; Donaldson, D. J.; Garrett, B. C.; George, C.; Griffiths, P. T.; Hanson, D. R.; et al. An Overview of Current Issues in the Uptake of Atmospheric Trace Gases by Aerosols and Clouds. Atmos. Chem. Phys. 2010, 10, 10561−10605. (20) Vaida, V. Perspective: Water Cluster Mediated Atmospheric Chemistry. J. Chem. Phys. 2011, 135, 020901. (21) Morokuma, K.; Muguruma, C. Ab Initio Molecular Orbital Study of the Mechanism of the Gas Phase Reaction SO3 + H2O: Importance of the Second Water Molecule. J. Am. Chem. Soc. 1994, 116, 10316−10317. (22) Vö hringer-Martinez, E.; Hansmann, B.; Hernandez, H.; Francisco, J. S.; Troe, J.; Abel, B. Water Catalysis of a RadicalMolecule Gas-Phase Reaction. Science 2007, 315, 497−501. (23) Lagutschenkov, A.; Fanourgakis, G. S.; Niedner-Schatteburg, G.; Xantheas, S. S. The Spectroscopic Signature of the “All-Surface” to “Internally Solvated” Structural Transition in Water Clusters in the n = 17−21 Size Regime. J. Chem. Phys. 2005, 122, 194310. (24) Lin, S.-S. Detection of Large Water Clusters by a Low rf Quadrupole Mass Filter. Rev. Sci. Instrum. 1973, 44, 516−517.

(25) Searcy, J. Q.; Fenn, J. B. Clustering of Water on Hydrated Protons in a Supersonic Free Jet Expansion. J. Chem. Phys. 1974, 61, 5282−5288. (26) James L. Kassner, J.; Hagen, D. E. Comment on “Clustering of Water on Hydrated Protons in a Supersonic Free Jet Expansion”. J. Chem. Phys. 1976, 64, 1860−1861. (27) Sloan, E. D. Clathrate Hydrates of Natural Gases; Marcel Dekker: New York, 1998. (28) Wei, S.; Shi, Z.; Castleman, A. W., Jr. Mixed Cluster Ions As a Structure Probe: Experimental Evidence for Clathrate Structure of (H2O)20H+ and (H2O)21H+. J. Chem. Phys. 1991, 94, 3268−3270. (29) Xantheas, S. S. Low-Lying Energy Isomers and Global Minima of Aqueous Nanoclusters: Structures and Spectroscopic Features of the Pentagonal Dodecahedron (H2O)20 and (H3O)+(H2O)20. Can. J. Chem. Eng. 2012, 90, 843−851. (30) McDonald, S.; Ojamäe, L.; Singer, S. J. Graph Theoretical Generation and Analysis of Hydrogen-Bonded Structures with Applications to the Neutral and Protonated Water Cube and Dodecahedral Clusters. J. Phys. Chem. A 1998, 102, 2824−2832. (31) Mortensen, J. J.; Hansen, L. B.; Jacobsen, K. W. Real-Space Grid Implementation of the Projector Augmented Wave Method. Phys. Rev. B 2005, 71, 035109. (32) Enkovaara, J.; Rostgaard, C.; Mortensen, J. J.; Chen, J.; Dułak, M.; Ferrighi, L.; Gavnholt, J.; Glinsvad, C.; Haikola, V.; Hansen, H. A.; et al. Electronic Structure Calculations with GPAW: A Real-Space Implementation of the Projector Augmented-Wave Method. J. Phys.: Condens. Matter 2010, 22, 253202. (33) TURBOMOLE V6.3 2011, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989−2007, TURBOMOLE GmbH, since 2007; available from http://www. turbomole.com. (34) Kuo, J.; Coe, J.; Singer, S.; Band, Y.; Ojamäe, L. On the Use of Graph Invariants for Efficiently Generating Hydrogen Bond Topologies and Predicting Physical Properties of Water Clusters and Ice. J. Chem. Phys. 2001, 114, 2527−2540. (35) Gregory, J. K.; Clary, D. C.; Liu, K.; Brown, M. G.; Saykally, R. J. TheWater Dipole Moment in Water Clusters. Science 1997, 275, 814− 817. (36) Anick, D. J. Application of Database Methods to the Prediction of B3LYP-Optimized Polyhedral Water Cluster Geometries and Electronic Energies. J. Chem. Phys. 2003, 119, 12442−12456. (37) Kirov, M. V.; Fanourgakis, G. S.; Xantheas, S. S. Identifying the Most Stable Networks in Polyhedral Water Clusters. Chem. Phys. Lett. 2008, 461, 180−188. (38) Eigen, M. Proton Transfer, Acid−Base Catalysis, and Enzymatic Hydrolysis. Part I: Elementary Processes. Angew. Chem., Int. Ed. Engl. 1964, 3, 1−19. (39) Marx, D. Proton Transfer 200 Years after von Grotthuss: Insights from Ab Initio Simulations. ChemPhysChem 2006, 7, 1848− 1870. (40) Marx, D.; Chandra, A.; Tuckerman, M. E. Aqueous Basic Solutions: Hydroxide Solvation, Structural Diffusion, and Comparison to the Hydrated Proton. Chem. Rev. 2010, 110, 2174−2216. (41) Marx, D.; Tuckerman, M. E.; Hutter, J.; Parrinello, M. The Nature of the Hydrated Excess Proton in Water. Nature 1999, 397, 601−604. (42) Hassanali, A. A.; Cuny, J.; Ceriotti, M.; Pickard, C. J.; Parrinello, M. The Fuzzy Quantum Proton in the Hydrogen Chloride Hydrates. J. Am. Chem. Soc. 2012, 134, 8557−8569. (43) Tuckerman, M. E.; Laasonen, K.; Sprik, M.; Parrinello, M. Ab Initio Simulations of Water and Water Ions. J. Phys.: Condens. Matter 1994, 6, A93. (44) Swanson, J. M. J.; Simons, J. Role of Charge Transfer in the Structure and Dynamics of the Hydrated Proton. J. Phys. Chem. B 2009, 113, 5149−5161. (45) Desiraju, G. R. Hydrogen Bridges in Crystal Engineering: Interactions without Borders. Acc. Chem. Res. 2002, 35, 565−573. (46) Santra, B.; Michaelides, A.; Fuchs, M.; Tkatchenko, A.; Filippi, C.; Scheffler, M. On the Accuracy of Density-Functional Theory 10834

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835

The Journal of Physical Chemistry A

Article

Exchange-Correlation Functionals for H Bonds in Small Water Clusters. II. The Water Hexamer and van der Waals Interactions. J. Chem. Phys. 2008, 129, 194111. (47) Kelkkanen, A. K.; Lundqvist, B. I.; Nørskov, J. K. Density Functional for van derWaals Forces Accounts for Hydrogen Bond in Benchmark Set of Water Hexamers. J. Chem. Phys. 2009, 131, 046102. (48) Wales, D. J.; Doye, J. P. K.; Dullweber, A.; Naumkin, F. Y.; Hodges, M. P. Global Minima of Protonated Water Clusters. http:// www-wales.ch.cam.ac.uk/CCD.html (2012). (49) Buch, V.; Milet, A.; Vácha, R.; Jungwirth, P.; Devlin, J. P. Water Surface Is Acidic. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 7342−7347. (50) Lee, K.; Murray, E. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy van der Waals Density Functional. Phys. Rev. B 2010, 82, 081101. (51) Wang, F.-F.; Jenness, G.; Al-Saidi, W. A.; Jordan, K. D. Assessment of the Performance of Common Density Functional Methods for Describing the Interaction Energies of (H2O)6 Clusters. J. Chem. Phys. 2010, 132, 134303. (52) Bahn, S. R.; Jacobsen, K. W. An Object-Oriented Scripting Interface to a Legacy Electronic Structure Code. Comput. Sci. Eng. 2002, 4, 56−66. (53) Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (54) Møller, C.; Plesset, M. S. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 618−622. (55) Weigend, F.; Ahlrichs, R. Balanced Basis Sets of Split Valence, Triple Zeta Valence and Quadruple Zeta Valence Quality for H to Rn: Design and Assessment of Accuracy. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (56) Rappoport, D.; Furche, F. Property-Optimized Gaussian Basis Sets for Molecular Response Calculations. J. Chem. Phys. 2010, 133, 134105. (57) Häser, M.; Ahlrichs, R. Improvements on the Direct SCF Method. J. Comput. Chem. 1989, 10, 104−111. (58) Hättig, C.; Weigend, F. CC2 Excitation Energy Calculations on Large Molecules Using the Resolution of the Identity Approximation. J. Chem. Phys. 2000, 113, 5154−5161. (59) Hättig, C.; Hellweg, A.; Kohn, A. Distributed Memory Parallel Implementation of Energies and Gradients for Second-Order Moller− Plesset Perturbation Theory with the Resolution-of-Theidentity Approximation. Phys. Chem. Chem. Phys. 2006, 8, 1159−1169. (60) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (61) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Atoms, Molecules, Solids, And Surfaces: Applications of the Generalized Gradient Approximation for Exchange and Correlation. Phys. Rev. B 1992, 46, 6671− 6687. (62) Hamann, D. R. H2O Hydrogen Bonding in Density-Functional Theory. Phys. Rev. B 1997, 55, R10157−R10160. (63) Xantheas, S. S. Ab Initio Studies of Cyclic Water Clusters (H2O)n, n = 1−6. III. Comparison of Density Functional with MP2 Results. J. Chem. Phys. 1995, 102, 4505−4517. (64) Lovas, F. J. Microwave Spectral Tables II. Triatomic Molecules. J. Phys. Chem. Ref. Data 1978, 7, 1445−1750.

10835

dx.doi.org/10.1021/jp307608k | J. Phys. Chem. A 2012, 116, 10826−10835