(H2O)20 Water Clusters at Finite Temperatures - American Chemical

Jun 3, 2013 - clusters and in the infrared spectrum at finite temperatures. In ... part, vibrational zero-point corrections, vibrational entropy, and ...
0 downloads 0 Views 691KB Size
Article pubs.acs.org/JPCA

(H2O)20 Water Clusters at Finite Temperatures 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: We have performed an exhaustive study of energetics of (H2O)20 clusters. Our goal is to study the role that various free-energy terms play in this popular model system and see their effects on the distribution of the (H2O)20 clusters and in the infrared spectrum at finite temperatures. In more detail, we have studied the electronic ground-state structure energy and its long-range correlation (dispersion) part, vibrational zero-point corrections, vibrational entropy, and proton configurational entropy. Our results indicate a delicate competition between the energy terms; polyhedral water clusters are destabilized by dispersion interaction, while vibrational terms (zero-point and entropic) together with proton disorder entropy favor them against compact structural motifs, such as the pentagonal edge- or face-sharing prisms. Apart from small water clusters, our results can be used to understand the influence of these energy terms in water/ice systems in general. We have also developed energy expressions as a function of both earlier proposed and novel hydrogen-bond connectivity parameters for prismatic water clusters.



for the magic number behavior at n = 21,3,11 but the role that the isomerism plays in the optical spectra has not been wellestablished, although some studies of this subject exist.38 In the present work, we take as a model system the neutral 20-molecule water cluster and study how vibrational and proton configurational entropic terms affect the distribution of structural motifs at finite temperatures. We study each of the four low-lying (H2O)20 families found in semiempirical potential studies39 and take into account all low-energy proton configurational isomers of each morphology with the aid of graph-theoretical techniques.33,37 This leads us to the concept of “proton configurational density of states”, a quantity that has a substantial contribution to the entropy of each structural motif and to the distribution of clusters. We also take into account the zero-point energies and vibrational entropy within the harmonic approximation. Two popular dispersion-corrected density functional theory (DFT) functionals are tested for water systems. We also study how the proton configurational density of states affects the infrared spectrum of (H2O)20. Another subject of this study is to derive energy expressions for (H2O)20 clusters using the so-called hydrogen-bond connectivity (HBC) parameters (a.k.a. topological invariants).33,37,40,41 With the aid of these parameters, cluster energies can be obtained simply from topological information of the hydrogen-bonded network without the need for explicit (say, semiempirical or ab initio) calculations. HBCs can be traced down to graph-invariants,37 that is, parameters that are

INTRODUCTION Small protonated water clusters ((H2O)nH+, n ≤ 21) have been studied for decades both experimentally1−3 and theoretically.4−14 Water clusters in general are relevant in many fields of science. In biology, small water clusters are important in intraprotein proton tranportation.15 In reactions occurring in the atmosphere, water is often present, as many of the reactions occur on small water droplets and aerosol surfaces.16−18 In chemistry, water in any quantity is known to promote many reactions by stabilizing the transition state.18−20 Especially the neutral (H2O)20 clusters have gained large interest4,21−31 because it has been speculated that the minimum-energy cluster would work as a building block for the “magic number” protonated cluster (H2O)21H+, which has been found to have accentuated stability in mass spectroscopy studies.1 Clusters of this size have also been used as a model of cloud droplets in atmospheric reactions.32 Neutral water clusters (H2O)n have been popular theoretical model systems for investigating fundamental aspects of hydrogen bonding and isomerism in water systems.4,25,30,32−35 Small hydrogen-bonded water clusters (n < 21) have several energetically low-lying structural motifs (polyhedral cages, prisms, cubes), while each motif can possess a large number of proton configurational isomers corresponding to rotational configurations of water monomers in the hydrogen-bonded network. Because of this, searching the proton configurational space is a challenging task for any potential energy surface (PES) sampling scheme, and therefore, a graph-theoretical technique has been developed to overcome this difficulty.33,36,37 Both types of isomerisms (morphological and proton configurational) should also be present in experiments involving protonated water clusters as the clusters retain finite internal energy. In fact, entropic effects have been speculated to account © 2013 American Chemical Society

Special Issue: Oka Festschrift: Celebrating 45 Years of Astrochemistry Received: January 10, 2013 Revised: May 31, 2013 Published: June 3, 2013 9985

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

electronic parts of the heat capacity at the constant volume (CV,trans, CV,rot, CV,vib, CV,elec) and the Boltzmann constant (k)

invariant under symmetry operations on the hydrogen-bond topology. In general terms, physical properties can be expressed as a function of the invariants. For example, expressing the heat capacity using graph-invariants makes it possible, in principle, to the study of phase transitions in ice systems.37 HBC parameters for polyhedral water clusters (PWCs) have been pioneered by Anick,30,34,35,40,41 while in the present work, we extend this approach to prismatic and other more “bulklike” clusters. We hope that in the future, this effort facilitates the screening of low-energy isomers for large and semi-infinite ice models and even the possibility of modeling the order− disorder phase transition of ice-Ih.

C P = k + C V,trans + C V,rot + C V,vib + C V,elec

The entropy at a temperature T and pressure P is calculated as a sum of translational, rotational, vibrational, and electronic entropies, that is S(T , P) = Strans + Srot + Selec + Svib − k × ln

Ni di exp( −Gi /kT ) = ∑j dj exp( −Gj /kT ) N

(1)

RESULTS Features of Cluster Morphologies. In this work, we consider various morphologies available for the 20-molecule water cluster, namely, the regular dodecahedron, edge-sharing pentagonal prisms (ESPPs), face-sharing pentagonal prisms (FSPPs), and fused cubes (FCs). These morphologies, which are depicted in Figure 1, have been established39 to be the most stable ones for a neutral water cluster of this size. Although these feature the same amount of water molecules, the morphologies differ from each other when considering the shape and the hydrogen-bonding topology: the dodecahedron is a PWC containing only three-fold-coordinated water molecules, while ESPPs and FSPPs consist of three pentagonal prisms that share an edge in ESPPs and a facet in FSPPs. The alignment of the prisms makes a difference when considering

(2)

T

C P dT

(6)



where H, T, and S are the enthalpy, temperature, and entropy of the system, respectively. The enthalpy and entropy are calculated using the ideal gas assumption with the IdealGasThermo module of ASE.42,50 In this model, the enthalpy at a given temperature T is calculated as

∫0

(5)

The previous equation (eq 6) is used to study the effect of proton configurational entropy on the distribution of water cluster morphologies. Different isomers and morphologies may be connected via high-energy pathways, which implies that the distribution is not necessarily ergodic. However, high-energy isomers are most likely produced during experimental procedures (supersonic jet expansion10,51−53), which are then quenched to local minima, that is, all local minima should be accessible and populated despite the high energetic barriers separating them. The synthesized clusters are also known to host considerable residual internal energy (35 kcal/mol),10 while temperatures of the clusters have been estimated to be ∼170 K.51

where E[X] is the electronic energy of species X. In order to study the energetic ordering of structures at finite temperatures, we evaluate the Gibbs free energies (G)

H(T ) = E + EZPE +

P(o)

In the previous equation, P is the pressure of the atmosphere (P(o) = 101.325 kPa). In this work, the pressure P is always equal to P(o). The vibrational frequencies needed to calculate the zero-point vibrational energies and vibrational parts of the heat capacity and entropy were computed using the harmonic approximation with the ASE module Vibrations.42 The infrared spectrum intensities were calculated “from a finite difference approximation of the gradient of the dipole moment” by the InfraRed module of ASE.42 Zero-point energy and vibrational terms were computed for a single (lowest-energy) isomer within each cluster morphology. For clusters in the proton configurational ensemble, the zero-point energy and entropy of eq 5 from the lowest-energy proton configurational isomer were used. The approximation was motivated by ref 35, which states that the zero-point energy is inversely related to the electronic energy so that a 1 kcal/mol increase in the electronic energy lowers the zero-point energy by 0.11 kcal/mol. The probabilities of finding a cluster i with a given free energy Gi value and a degeneracy factor di are obtained from the Maxwell−Boltzmann distribution law

METHODS All structure relaxations were performed using ASE (Atomic Simulation Environment) version 3.4.42 Two different codes were used for electronic structure calculations; DFT calculations were performed with GPAW version 0.843,44 and MP2 (second order Møller-Plesset perturbation theory) calculations with the TurboMole program package.45 Structural relaxations for large structural samples were terminated when the residual force was lower than 0.05 eV/Å, while for single structures, a tolerance of 0.02 eV/Å was used (1 eV = 96.485 kJ/mol, 1 Å = 1 × 10−10 m). In GPAW calculations, the grid spacing was 0.2 Å, and the systems were relaxed in a supercell having 5 Å of surrounding vacuum in each direction. The TurboMole MP2 calculations were performed with a def2-TZVPD (triple-ζ valence with polarization and diffuse functions) basis set. In the DFT calculations, three different exchange and correlaction (XC) functionals were used, namely, the PBE,46 vdw-DF2,47 and TS09 functionals.48 The latter two include dispersion corrections; vdw-DF2 relies on a two-body longrange correlation kernel, evaluated either in the real or Fourier space, while TS09 uses Grimme’s dispersion correction scheme,49 with “on-the-fly” corrected C6 terms. Due to the cost of evaluating the long-range correlation term in vdw-DF2, the TS09 correction is computationally more efficient. The conventional parts of the XC correlational functionals in our calculations were PBE and PW86 for TS09 and vdW-DF2, respectively (following refs 47 and 48). Stabilization energies ΔE of neutral clusters (H2O)n are calculated as follows

G (T ) = H (T ) − T × S(T , P )

P

(o)



ΔE = E[(H 2O)n ] − nE[H 2O]

(4)

(3)

where EZPE is the zero-point vibrational energy. The quantity CP is the heat capacity at constant pressure. It can be calculated as a sum of the translational, rotational, vibrational, and 9986

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

is always ∼108°, which is well-suited for achieving the linear O−H−O bond, but on the other hand, its constituent water monomers are only three-fold-coordinated and thus less hydrogen-bonded. The FC and FSPP morphologies can also be classified as INTs,28 which form vertical molecular rings placed on top of each other; an FC can be classified as an INT that contains five layers, each containing four water molecules, while an FSPP is classified as an INT structure with four layers of five molecules. In this work, we also study some smaller INTs, denoted as INTmn , where m is the number of layers and n the number of molecules in a layer. The structures to be studied are INT44, INT43, and INT35 (Figure 1). The motivation for working on these smaller INTs is to investigate if the energetics of larger INTs can be characterized by HBC parameters that were fitted using smaller INTs. Our objective is to study the energetic stability of (H2O)20 morphologies at finite temperatures, once the proton disorder, dispersion correction, zero-point energies, and vibrational entropy are taken into account. To our knowledge, there is no earlier comprehensive study of (H2O)20 that would take into account all of these energy contributions. There is a widespread consensus that ESPP is the most stable neutral 20-molecule water cluster at 0 K and likely at cryogenic temperatures as well,21−26 although in some studies, FSPP has been found to be to the most stable species.28,29 However, the assumption that ESPP is the prominent species at cryogenic temperatures is based on ab initio calculations of single-proton configurations of particle morphologies, omitting proton configurational entropy, and, in some cases, zero-point energies, which are also important for 0 K energy ordering and for other freeenergy factors. A study by Lenz et al. states that the proton configurational entropy makes the dodecahedron the most common species above 200 K, but in that paper, the conclusions are made with many approximations, and FSPP and FC morphologies are not taken into account.27 Structure Generation. All hydrogen-bond topologies for each (H2O)20 and INT species were generated using our own program package that exploits the graph invariants introduced by Kuo et al.37 in finding the symmetry-independent clusters. The number of topologies found is tabulated for each morphology in Table 2. The numbers agree with previous findings.28,33,37,54 The largest number of topologies is found for the dodecahedron as it has the smallest number of 4-fold coordinated molecules (i.e. C-type molecules, see Table 5) molecules, while most of the molecules in the FC cluster are four-fold-coordinated, yielding the smallest number of top-

Figure 1. Most stable hydrogen-bond topologies of (H2O)20’s and smaller INTs. Black lines are hydrogen bonds, and red and white circles denote oxygen and hydrogen atoms, respectively.

the hydrogen-bonding pattern, as is evident from Table 1. The oxygen atom “raft”, which is unique in all morphologies, creates Table 1. Features of Studied Cluster Morphologiesa topology

NC

NF

NL

symmetry group

max αO (deg)

min αO (deg)

dodecahedron ESPP FSPP (INT45) FC (INT54) INT35 INT44 INT34

0 8 10 12 5 8 4

10 6 5 4 5 4 4

10 6 5 4 5 4 4

Ih D3h D5h D4h D5h D4h D4h

108 160 180 180 180 180 180

108 90 90 90 90 90 90

a Different monomer species are defined in Table 5. See Table 5 for the meanings of used abbreviations. The angles are given in degrees.

hydrogen bonds of different quality. In ice nanotubes (INTs), the optimal O−H−O angle for the hydrogen bonding, that is, 180°, can be hard to achieve as there are O−O−O angles αO of ∼90 and ∼180°. The ESPP also has some αO angles that are unfavorable, namely, ∼90 and ∼120°. In the dodecahedron, αO Table 2. Features of the Studied (H2O)20 Clustersa topology

ntotal

nreduced

nsampled

min dO−O (Å)

max dO−O (Å)

avg dO−O (Å)

min ΔE avg dO−O (Å)

dodecahedron ESPP FSPP (INT45) FC (INT54) INT35 INT44 INT34

30026 22960 7480 5588 860 978 178

30026 4537 408 255 204 134 66

1458 1282 408 255 204 134 178

2.70 2.76 2.77 2.77 2.77 2.76 2.77

2.85 2.83 2.84 2.85 2.84 2.85 2.84

2.76 2.79 2.79 2.80 2.80 2.80 2.81

2.72 2.76 2.77 2.78 2.77 2.76 2.81

a

The quantity ntotal denotes the total number of hydrogen-bond topologies, nreduced the number of topologies with no unfavorable O−O−O angles (see text), and nsampled the number of topologies optimized with PBE+TS09. Additionally, oxygen distances (Å) between nearest-neighbor water molecules of neutral ice clusters are tabulated. Min dO−O, max dO−O, and avg dO−O denote the minimum, maximum, and average O−O distance, respectively. Min ΔE avg dO−O is the average O−O distance of the minimum-energy structure. 9987

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

Figure 2. The stabilization energy distributions (i.e., proton configurational density of states) of (H2O)20 morphologies (the dodecahedron [blue], ESPP [green], FSPP [red], and FC [teal]) calculated with (a) PBE and (b) TS09 including nsampled structures (see Table 2).

Table 3. Minimum Stabilization Energies (per molecule in parentheses) of Different Morphologies of (H2O)20 Compared to ESPP Calculated with Different Levels of Theorya method

dodecahedron

MP2/def2-TZVPD PBE vdW-DF2 TS09

11.71 1.61 7.67 6.97

(0.59) (0.08) (0.38) (0.35)

MP2/aug-cc-pVTZ TS09

11.04 (0.56) 7.32 (0.37)

ESPP 0 0 0 0

[−220.82] [−211.38] [−206.70] [−233.56]

0 [−220.32] 0.52 (0.03)

FSPP 1.24 3.09 1.53 2.50

(0.06) (0.15) (0.08) (0.13)

1.87 (0.09) 2.89 (0.14)

FC 1.84 7.36 2.51 5.32

(0.09) (0.37) (0.13) (0.27)

2.44 (0.12) 6.04 (0.30)

a

The absolute stabilization energy of ESPP is marked in square brackets. Previously reported MP2 results26 together with TS09 values are given in the lowermost part of the table. All numbers are in kcal/mol.

of the dodecahedron and ESPP. Finally, all results were recalculated with the TS09 XC functional to study the contribution of the dispersion interaction. During most structural optimization runs, the structures remained close to the initial geometries. In some cases, considerable changes were observed for the dodecahedron: (1) of the 1458 dodecahedron topologies, in 76 structures, a water molecule spontaneously ionized, causing the formation of a H3O+ and OH− ion pair; (2) in 65 structures, an intramolecular OH bond was elongated so that one hydrogen atom was nearly in the midpoint between two oxygen atoms, that is, forming a Zundel-like species. Kuo et al. made a similar observation regarding the spontaneous ionization of some cluster geometries in ref 55. Average O−O distances between nearest-neighbor water molecules for different structures are tabulated in Table 2. The values are similar between the structures. The only small variation occurs for the dodecahedron, where the bond lengths are slightly shorter (0.02−0.06 Å) than those for the other structures. It is noteworthy that the numbers are almost equal for all INTs, and no variance of the O−O bond length as a function of the number of rings in the structure is observed, which is in contrast to the observations reported in ref 28. Electronic Structure Energetics. Within the framework of DFT, all studied cluster structures were calculated using both the PBE and TS09 XC functionals. The effect that vdW corrections have on the proton configurational density of states is clear in Figure 2; dodecahedral structures are shifted up in stabilization energy as dispersion favors compact structures. The effect of dispersion was further tested with several computational methods for the most favorable isomers as

ologies. The number of hydrogen-bond topologies for the dodecahedron would be even larger if it did not possess such a high symmetry. Among the automatically generated topologies, that is, when hydrogen atoms are placed into the oxygen atom raft, there are many clusters containing unnaturally large intramolecular H− O−H bond angles (over 120°, as discussed in the previous section). If these unnatural cluster structures are dropped out, the number of different hydrogen-bond topologies is significantly reduced for most of the clusters to nreduced structures (tabulated in Table 2). We started our study of cluster energetics by relaxing some of the topologies using DFT calculations and the PBE exchange and correlation functional (performed for ESPP in our previous study4). In the calculations, we considered only the physically meaningful nreduced clusters. However, to check the validity of using only nreduced number of clusters, all ntotal structures were calculated for INT34. Unnatural bond angles resulted in weaker hydrogen bonds and less favorable stabilization energies, confirming that dropping them out is a reasonable approximation. For the morphologies containing (tens of) thousands of hydrogen-bond topologies (the dodecahedron and ESPP), all nreduced structures were not calculated, but a smaller set of nsampled structures was chosen. For the dodecahedron, we calculated nsampled = 1458 topologies, of which most were picked up randomly. Approximately 100 low-energy topologies, with the smallest possible bFF value (see Table 5), were included in this number. The same scheme was also used for ESPP. This procedure ensures a correct description of the proton configurational density of states at the energy minimum 9988

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

Figure 3. Relative stabilization energies of different (H2O)20 morphologies calculated using different levels of theory. The stabilization energy of ESPP with each method is set to zero. The symbol D indicates the dodecahedron.

Table 4. Relative Gibbs Free Energies (G), Divided into the Entalphy (H) and Entropy (−TS) Contributions for the MinimumEnergy Hydrogen-Bond Topologies of Table 3a dodecahedron

ESPP

FSPP

FC

T/K

H

−TS

G

H

−TS

G

H

−TS

G

H

−TS

G

273.15 253.15 233.15 213.15 193.15 100 0

3.87 3.83 3.78 3.74 3.69 3.44 3.18

−2.05 −1.87 −1.69 −1.51 −1.32 −0.49 0.00

1.80 1.94 2.10 2.24 2.38 2.93 3.18

0 (19.12) 0 (16.44) 0 (13.91) 0 (11.55) 0 (9.39) 0 (2.10) 0

0 (−49.15) 0 (−42.96) 0 (−37.14) 0 (−31.70) 0 (−26.66) 0 (−8.85) 0

0 (−30.02) 0 (−26.52) 0 (−23.22) 0 (−20.15) 0 (−17.27) 0 (−6.76) 0

1.66 1.55 1.50 1.41 1.31 0.95 0.83

−1.45 −1.28 −1.11 −0.95 −0.78 −0.14 0.00

0.18 0.28 0.37 0.46 0.53 0.81 0.83

5.88 5.86 5.86 5.86 5.86 6.20 6.60

2.27 2.11 1.95 1.78 1.59 0.56 0.00

8.14 7.98 7.79 7.63 7.45 6.78 6.60

a

Energies are expressed in the units of kcal/mol, and at each temperature, numbers relative to the corresponding ESPP values are listed. Absolute values for the ESPP case are in parentheses. The electronic structure energy (E) and zero-point energy (EZPE) of eq 3 have been set to zero for ESPP (EZPE = 326.71 kcal/mol); the row corresponding to T = 0 K shows the relative energy spacing of structures with zero-point corrections included (values that can be compared to the results in Table 3).

for MP2. This could be due to the well-known overestimation of dispersion in MP2.56,57 Energy spacings among the compact ESPP, FSPP, and FC structures show a spread between the different dispersioncorrected functionals, especially for the FC case, which is more stable (and closer to MP2 results) with vdW-DF2 than that with TS09. The agreement between the dodecahedron, ESPP, and FSPP is reasonable among the dispersion-corrected DFT functionals (≤1 kcal/mol). Calculations with vdW-DF, MP2, and CCSD(T) (the coupled cluster method including single and double excitations fully and triple excitations noniteratively) have been compared in ref 58, and a spread of ∼0.1 kcal/mol per molecule was observed among different methods for a certain water hexamer isomer. In the present case (20 water molecules), we may estimate a roughly ≤2 kcal/mol error in the relative sublimation energies. Finite Temperature. In this section, we first consider the effect of the translational, rotational, and vibrational entropy on the energetics of 20-molecule water clusters. After this, we take into account the proton disorder entropy, which originates from the proton configurational isomers of each water cluster species. The vibrational frequencies for minimum-energy structures were calculated at the harmonic approximation and by exploiting the TS09 XC functional. These were used to compute the free energy as a function of the temperature with

identified using the TS09 XC functional. Results of this test are listed in Table 3 and visualized in Figure 3. From the minimum stabilization energies, it is observed that our MP2 results are similar to the previous ones as given by Fanourgakis et al.26 (also shown in Table 3). The energies for ESPP, FSPP, and FC are only 0.50, 0.74, and 1.10 kcal/mol smaller than those in ref 26, and the energy ordering (ESPP < FSPP < FC < the dodecahedron) of the morphologies is unchanged. The smaller energy spacing is caused from use of distinct hydrogen-bond topologies compared with ref 26 (those topologies were taken from the semiempirical potential calculations of Wales and Hodges39). According to the results of Table 3, our minimum-energy structures, as found with a full graph-theoretical search, are slightly more optimal (by 0.50 kcal/mol for ESPP) than those in ref 26. Comparing relative stabilization energies in Table 3 and Figure 3 shows clearly how important the dispersion corrections are in order to describe the energetic ordering between “hollow” (the dodecahedron) and “compact” (prismatic) structures correctly. For example, the difference between ESPP and the dodecahedron is 1.61 kcal/mol at the PBE level, while with MP2, it is 11.71 kcal/mol. Dispersioncorrected DFT functionals (vdW-DF2 and TS09) are, according to Table 3 and Figure 3, in qualitative agreement with MP2 in terms of energy spacing but give a slighty smaller energy spacing between ESPP and the dodecahedron than that 9989

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

Figure 4. Fractions of (H2O)20 morphologies calculated with eq 6 at different temperatures using the (a) electronic energy, (b) free energy, and (c) free energy with the proton configurational density of states (Figure 5). The dodecahedron, ESPP, and FSPP have been marked with blue, red, and green colors, respectively.

Figure 5. Distribution of proton configurations near the minimum energy at different temperatures. ESPP, FSPP, and the dodecahedron are marked with green, red, and blue, respectively (the free energy of the minimum-energy hydrogen-bond topology of ESPP is set to zero).

is included. A similar change in relative energies due to ZPEs has been previously found for the prism and cage morphologies of the water hexamer.59,60 The effect that the zero-point energies and entropy have on the distribution of cluster morphologies is depicted in Figure 4; panel (a) shows a case where only the electronic structure energy is considered, while panel (b) gives the distribution with the zero-point energy and entropy taken into account. Next, we will consider the effect that the proton configurational density of states of Figure 2 has on the distribution of cluster morphologies. Figure 2b depicts the configurational density of states when only the electronic structure energy is considered. Figure 5a shows how the distributions are shifted when the zero-point energy is taken into account. We notice that the minimum-energy topologies of FSPP move closer to

eq 2. Results are listed in Table 4, where relative contributions (with respect to ESPP) of each energy term have been tabulated. According to Table 4, even when the zero-point corrections and entropy are taken into account, ESPP is energetically the most stable cluster geometry at all temperatures. However, comparing the values of the T = 0 K row of Table 4 with the values of Table 3 for ESPP and FSPP shows that zero-point energies close up the energetic gap between these two morphologies from 2.50 to 0.83 kcal/mol. Including the temperature and taking into account the entropy closes this gap further, that is, at T = 193.15 K, it is only 0.53 kcal/mol. A similar effect is observed for the dodecahedron as taking into account the zero-point energy lowers its energy difference with ESPP from 6.97 to 3.18 kcal/mol and even more as the entropy 9990

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

of Figure 6 are similar as only slight changes in the shapes of the peaks are observed. This means that the significance of proton disorder is quite small for ESPP as it has very few other proton configurations with significant weights (see Figure 5). When all of the (H2O)20 proton configurations are taken into account (Figure 6c), large changes in the spectrum can be observed; new broader peaks of smaller intensity emerge near 2300, 2600, 2700, and 3100 cm−1. In addition, some intensity appears between 3250 and 3400 cm−1, the relative intensities of the peaks are altered, and some new sharp features appear. The broad features are caused by the tens of dodecahedron configurations that are only slighty different from each other, causing the appearance of many partially overlapping peaks. The few significant FSPP configurations cause the new sharper peaks at approximately 3150 and 3200 cm−1 and the intensity gain of the 3500 cm−1 peak. Energy Expressions with Connectivity Parameters. In this section, we discuss how the stabilization energies of clusters depend on the HBC parameters. The abbreviations used are explained in the section Hydrogen-Bond Connectivity Parameters in Appendix A. We concentrate on ESPP, FSPP, FC, and INTs as PWCs have been studied in detail earlier.30,33−35,40,41,55,61,62 Unnatural bond angles and structures that form ions spontaneously (see the section Structure Generation) are excluded. The abbreviations used for HBCs are denoted in Table 5 in Appendix A. Earlier investigations on PWCs gave a starting point for our work as each of the studied morphologies contains many L- and F-type molecules in addition to C-type molecules. As emphasized in previous studies,30 bFF is a significant factor in the electronic energy of PWCs. We studied if this is the case for clusters considered in the present work. For this purpose, we plotted the stabilization energy of the (H2O)20 clusters as a function of the dipole moment and grouped the results with the number of FF bonds for each morphology (Figure 7). From Figure 7, it is obvious that bFF plays a significant role in the electronic energy of the cubic and pentagonal morphologies, but it is not as remarkable as that for the dodecahedron. There is an energetic overlap of clusters having the same bFF values, especially for FSPP and FC. There are other important factors affecting the energetics as there exist several energetically separated clusters with the same bFF values in the figure. To find out these factors, several connectivity parameters were tried. We observed that an important parameter is bFC, similarly as that in (H2O)21, as studied in ref 30, but for morphologies studied in this article, the parameter is even more critical. Another significant factor turned out to be the number of threemolecule chains, that is, aTh. Grouping the results with bFF, bFC, and aTh causes the energy−dipole moment plot to form clear, energetically separated groups. This is visualized for FSPP and FC in Figure 8. An interesting finding is such that in Figure 8, in many cases, the energy is independent of the dipole moment. For example, an FSPP cluster can have a dipole moment of 6 or 15 D with practically the same energy, which is in contrasts to previous findings for the dodecahedron and ESPP.4,30,62 The smaller INTs also form clean groups with the same parameters, but for these morphologies, a significant energy−dipole moment dependence is observed. We further tuned the model described above with novel parameters (see Table 5). Of these, the parameter that was found to have the most significant effect on the energetics was aFFC. In addition, we used parameters si, pi, g4, and o5, which can

the minimum-energy topologies of ESPP. As the temperature and entropy are included, the shift grows larger, which is evident from Figure 5b,c. The distribution of clusters at various temperatures is depicted in Figure 4c. It is observed that even at temperatures as low as ∼200 K, a considerable amount of FSPP and dodecahedral morphologies is present. It must be emphasized that the distribution of cluster morphologies depends heavily on the accuracy of the used electronic structure calculation method, especially on the quality of the computed dispersion interaction. Looking at Figure 5, it is obvious that a downward shift of ∼1 kcal/mol in the energies of the dodecahedrons would saturate the lowenergy region of the proton configurational density of states completely, making the dodecahedron the most prominent cluster morphology even at small temperatures as it has the highest number of isomers. Also, as discussed in the previous sections, the relative positions of the dodecahedrons in Figure 5 depend on the amount of dispersion, which should be described with high accuracy. According to Table 3, the energy spacing between ESPP and FSPP may jump when changing from TS09 XC to MP2, ∼1 kcal/mol, which would, according to Figure 5, diminish the portion of ESPP at temperatures even at 100 K. The small energetic shifts among the morphologies also depend sensitively on the zero-point energy and entropy, both calculated in the present work using the dispersioncorrected DFT and harmonic approximation. We stress that although our results are qualitative, a large fraction of other morphologies apart from ESPP is present in the distribution. To illustrate the effect that the proton configurational density of states has in experimental measurements, we calculated benchmarking harmonic infrared absorption spectra of (H2O)20 in a few different cases at 233.15 K (Figure 6). In Figure 6a,

Figure 6. Infrared spectra of (H2O)20 at 233.15 K depicted (a) for the most stable (ESPP) proton configuration, (b) for all ESPP clusters, and (c) for all (H2O)20 proton configurations. In (b) and (c), the spectrum of each proton configuration is multiplied by Boltzmann weights calculated with eq 6 using the free energy. Intensities on the vertical axis are in arbitrary units.

only the spectrum of the most stable ESPP proton configuration is present. Figure 6b shows the spectrum where the contributions of other ESPP particles are also added with the weights obtained from the Boltzmann equation (eq 6) using free energies. The entire (H2O)20 proton configurational density of states is considered in Figure 6c. Panels (a) and (b) 9991

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

Figure 7. The stabilization energies of different (H2O)20 morphologies ((a) the dodecahedron, (b) ESPP, (c) FSPP, and (d) FC) plotted as a function of the dipole moment. The points are colored by the number of FF bonds (bFF) in the structure (1 D = 3.33564 × 10−30 C m).

minor. However, these give, as mentioned earlier, the final energy ordering inside of the colored groups of Figure 8. Dodecahedron:

be employed to fine-tune the energetics inside of the colored groups of Figure 8. In order to study if the total energy depends on the total dipole moment of the cluster, the square of the dipole moment (|M|2) of the cluster was also included in the fit. Only independent fitting parameters were used (dependencies among parameters are listed in Appendix A). For ESPP, FSPP, and INT35, aFFF and aLLL instead of aTh were used as the former parameters are, in this case, independent. Finally, we selected the parameter set so that the regression would give positive factors for most of the parameters. The regression formulas with Pearson’s R parameter and root-mean square error (RMSE) values found for each morphology are listed in eqs 7−13. The stabilization energies (ΔE) given by these formulas are in kcal/mol. From the formulas, several conclusions can be drawn: (1) The linear regression of the HBC parameters to the stabilization energies works well, as can be observed from the R and RMSE values (sub-kcal/mol accuracy). (2) The stabilization energies of the FSPP and FC topologies are practically independent of the dipole moment, while for the small INTs and ESPP, there is some energy− dipole moment correlation. (3) The significance of si and pi is

ΔE = 0.561836(bFF − 3) + 1.125428a Th + 0.116999p3 + 0.067147p2 + 0.071813p1 − 0.147702g4 + 0.049819o5 + 0.016901|M |2 − 226.980086 R = 0.99826,

RMSE = 0.334

(7)

ESPP: ΔE = 1.910316bFF + 0.885661bFC + 0.784390aFFF + 0.387230aLLL + 0.280186aFFC + 0.038667g4 − 0.449527s1 − 0.185070s2 − 0.174272s3 − 0.078066p1 − 0.063171p2 − 0.048869p3 + 0.011747|M |2 − 231.023290

R = 0.99666,

RMSE = 0.302 (8)

9992

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

Figure 8. The stabilization energies of (a) FSPP and (b) FC morphologies plotted as a function of the dipole moment. The points in the graph are colored by bFF, bFC, and aTh.

FSPP (INT54):

INT34: ΔE = 1.447571bFF + 1.058662bFC + 0.351081a Th

ΔE = 0.947569(bFF − 1) + 1.313831bFC + 1.486710aFFF + 1.332040aLLL + 0.463444aFFC + 0.585143o5

+ 0.328617aFFC + 0.053998s2 + 0.088350s3

+ 0.115753g4 + 0.058831s2 + 0.115205s3 + 0.909861p1

+ 0.214579g4 + 0.011862|M |2 − 129.643300

+ 1.016452p2 + 0.728801p3 + 0.001247|M |2

R = 0.99837,

− 231.660179R = 0.99668,

RMSE = 0.238

CONCLUSIONS Understanding the role of free energy in the distribution of hydrogen-bonded water clusters is important when interpreting experimental infrared (IR) spectra as various isomers can be present in experiments involving protonated water clusters at finite temperatures. However, in theoretical studies, typically spectra from a single isomer are only considered and compared to the experimental results.10 The role that various isomers play in the IR spectra of water clusters has not been wellestablished,2 while some experimental3 and theoretical11 works suggest that the magic number (H2O)nH+, n = 21, results from entropic effects. In this study, we present the IR spectrum of (H2O)20 where all energetically relevant proton configurations of all four morphologies are taken into account. We studied the effect of the zero-point energy, vibrational entropy, and proton configurational entropy in the distribution of a popular model system, the neutral 20-water molecule cluster. According to our calculations, while the edge-sharing prism motif (the most stable isomer when considering electronic structure energies only) is the dominant species at temperatures below 200 K, a considerable amount of other lowenergy motifs, namely, the dodecahedron and the face-sharing prism, are present. At T ≈ 200 K, they constitute approximatively half of the clusters under equilibrium conditions. As a result, the IR spectrum of (H2O)20 is very different compared to the spectrum of the most stable ESPP morphology. Both the zero-point energies and the vibrational entropy stabilize the face-sharing prism and dodecahedron motifs against the edge-sharing prism, while the dodecahedral morphology is destabilized due to van der Waals forces, the net

ΔE = 1.809284(bFF − 1) + 1.207851bFC + 0.277554aFFF + 0.225301aLLL + 0.388375aFFC + 0.073747o5 + 0.039335s2 + 0.073342s3 + 0.676866p1 + 0.620328p2 + 0.460871p3 + 0.013373|M |2 − 167.497156 RMSE = 0.155

(10)

FC (INT54): ΔE = 1.646736bFF + 1.310324bFC + 1.135847a Th + 0.512990aFFC + 0.084338s2 + 0.135469s3 − 0.128117g4 + 0.001148|M |2 − 228.417582 R = 0.99799,

RMSE = 0.197

(11)

INT44: ΔE = 1.404319bFF + 1.284482bFC + 0.246728a Th + 0.607133aFFC + 0.104895s2 + 0.154469s3 − 0.374169g4 + 0.011375|M |2 − 181.079873 R = 0.99644,

RMSE = 0.232

(13)



(9)

INT35:

R = 0.99787,

RMSE = 0.139

(12) 9993

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

effect being that it lies ∼2 kcal/mol above the most stable prismic isomer. However, as it has the highest proton configurational density of states near the energy minimum and at elevated temperatures, it quickly gains the largest entropic weight. Because of this, the IR spectrum peaks caused by the dodecahedron are very broad when compared to the ones resulting from vibrations of ESPP and FSPP particles. Considering all proton configurational isomers near the energy minimum reveals that the minimum becomes mixed with various proton configurational isomers of both the edgeand face-sharing pentagonal prism motifs within sub kcal/mol energy spacing. Fine-tuning the cluster distribution becomes tedious; ∼2 kcal/mol uncertainties are present in the MP2 calculations of 20-water monomer water clusters, and energetic shifts of this order can tip the distribution in favor of the facesharing prism or the dodecahedron as the relative energy positions are sensitive to the quality of the used approximation for dispersion, an issue that could only be resolved with the “gold-standard” electronic structure method, that is, with the computationally expensive coupled cluster calculations at the CCSD(T) level. Concluding, proportions of low-energy structural motifs due to zero-point energy and entropic terms, including proton disorder, were estimated at typical experimental temperatures. However, the fine details of the distribution depend sensitively on the employed computational scheme and its accuracy, especially on the dispersion interaction. This presents a challenge for global minima search schemes that are based on semiempirical energy calculations.8,11,39 The second major effort in this work was to develop expressions for electronic structure energies of small water clusters using the so-called hydrogen-bond connectivity parameters. Earlier, this was done in the context of polyhedral water clusters, while in the present work, we have extended this idea to cubic and prismatic clusters. Cluster energies were parametrized with connectivity parameters introduced earlier in the literature,30,33−35,37,40,41 namely, as a function of singlemolecule types and two- and three-molecule hydrogen-bonded chains, but also using novel parameters, such as classified square- and pentagon-type hydrogen-bonded loops. Using these parameters, energy fits for the proton configurational ensemble with sub-kcal/mol errors only were obtained. We hope that our results and further developments will, in the future, be used in screening low-energy isomers in larger hydrogen-bonded systems. Efficient parametrization of the energy and heat capacity of hydrogen-bonded systems also creates possibilities in studying order−disorder phase transitions in model water systems and in ice-Ih.



Table 5. Nomenclature of HBC Parameters C L F NX bXY aXYZ

aTh si pi g4 o5

water molecule that accepts and donates two hydrogen bonds (AADD); see Figure 9a water molecule that accepts one hydrogen bond and donates two (ADD); see Figure 9b water molecule that accepts two hydrogen bonds and donates one (AAD); see Figure 9c Number of molecules of type X number of two-molecule hydrogen-bonded chains (pairs) where the donor is of type X and the acceptor is of type Y (see Figure 9d,c) number of three-molecule hydrogen-bonded chains, where the first, second, and third molecules are of type X, Y, and Z, respectively; we consider only aFFF, aLLL, and aFFC chains (the aFFC case is not unique; we consider only the case where the middle F molecule accepts a hydrogen bond from F and donates one to the C molecule total number of single-type three-molecule (i.e., FFF or LLL) chains, i.e., aTh = aFFF + aLLL number of hydrogen-bonded square loops of type i (types are defined in Figure 10) number of hydrogen-bonded pentagon loops of type i (types are defined in Figure 11) number of hydrogen-bonded square or pentagon loops with either four F or L molecules Number of hydrogen-bonded square or pentagon loops with either five F or L molecules

studies,30,34,35,40,41,61 while aFFC, pi, and sj have been introduced by us. In the following, we give a brief rationale for using the parameters discussed above. Suppose we have a term bXYk in the expression for the energy, where bXY is the number of twomolecule chains (see Table 5) and k is the parameter to be optimized. We rewrite this term as follows b XY

bXY k =

∑k = ∑ i=1

(i , j)XY

C(i , j) (14)

On the-right hand side, the summation indices i and j designate different molecules in the structure, the summation (i,j)XY goes over neighboring molecule pairs that are of type X and Y, and C(i,j) describes the interaction between such molecules. The trivial sum term in the middle running from 1 to k is given to emphasize that there are bXY times the k interaction term and to make a connection between the leftmost and rightmost expressions. The same model can be extended to N-molecule hydrogen-bonded chains. One then sums up the interaction energies between monomers, not explicitly between all of them but along hydrogen-bonded chains. This can be contrasted to the standard procedure (say, TIP4/5P water models) to calculate the electrostatic interaction energy, where the rightmost sum in eq 14 would run over all possible atomic pairs, that is, not only molecular pairs connected by hydrogen bonds. The fitting of ab initio results into an equation of type eq 14 has advantages as the fitting takes implicitly into account charge-transfer effects, with a minimal number of terms calculated. In the present work, we also take into account the energy dependence (and implicitly, the charge transfer) in some typical structural motifs (squares and pentagons) such as those depicted in Figures 10 and 11. Once the fitting parameters are known, energetics of proton disorder can be accessed at minimal computational cost. One merely has to generate the structures, while no total energy calculations and/or structural relaxations are needed. Using small model systems, such as the clusters considered in this study, parameters can be optimized using ab initio energy values as data, and after, this they can, in principle, be recycled

APPENDIX A

Hydrogen-Bond Connectivity Parameters

One of the objectives of this work is to derive expressions for the electronic energy of hydrogen-bonded systems as a function of the so-called HBCs. These have been previously developed for PWCs by Anick30,34,35,40,41 and Lenz et al.,61 while the error in PWC energy when evaluated by this method for a small set of isomers is only a fraction of a kcal/mol (1 kcal/mol = 4.184 kJ/mol).30 The HBC parameters used in this work are defined in Table 5 and Figure 9. The nomenclature is similar to the one introduced in ref 40. Some equations relating different HBCs for the cluster morphologies studied in this work are given. Most of the parameters have been used in earlier 9994

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

Figure 9. Water molecules with different hydrogen-bond acceptor/donor character: (a) double-acceptor, double-donor “C”-type, (b) singleacceptor, double-donor “L”-type, and (c) double-acceptor, single-donor “F”-type species. FL- and LC-type hydrogen bonds are depicted in panels (d) and (e), respectively.

Figure 10. Available hydrogen-bonding patterns in a four-molecule hydrogen-bonded loop (si). The red and black circles denote oxygen and hydrogen atoms, respectively. These loops are found in all studied clusters except in the dodecahedron.

Figure 11. Available hydrogen-bonding patterns in a five-molecule hydrogen-bonded loop (pi). The red and black circles denote oxygen and hydrogen atoms, respectively. These loops are found in the dodecahedron, ESPP, FSPP, and INT35 (see Figure 1).

in calculating energetics of larger systems. This approach can be powerful for semi-infinite systems (slabs) as the structures become rapidly too large for ab initio calculations when increasing the cell size. In the present work and in the context of HBCs, we take a step from PWCs to the direction of prismatic and more compact and bulk-like structures. Once energy expressions are available for this class of systems, lowenergy candidates can be efficiently screened for performing ab initio calculations. Finally, we emphasize that HBCs are invariant in the sense that they are not altered under symmetry operations on the hydrogen-bond network. Expressions between the “intuitive” parameters of Table 5 and “formal” graph-invariants can be derived. For example, Kuo et. al. have derived this kind of expression for NX for the (H2O)6 cluster.37

Relations between HBC Parameters

In this section, relations between connectivity parameters (see Table 5) are derived in order to find independent parameters in ice-rule-obeying INTs and prismic clusters. Dependencies between parameters have previously been derived for PWCs (including the dodecahedron) by Anick.40,41 The first connectivity parameter to be considered is obviously the number of hydrogen bonds in the cluster (Nb). As the number of C-type molecules is different in each water cluster species, the number of hydrogen bonds varies. The total number of bonds Nb for each species is the total number of hydrogens in the cluster (twice the number of water molecules NH2O) minus the number of dangling hydrogen bonds, which in our case is 9995

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

the number of F-type molecules (NF). Then, the following holds: Nb = 2 × NH2O − NF

s4 = s1 + 5

s4 = −0.5(s2 + s3) + 13

4

∑ si = 4(m − 1) + m

(15)

(24)

i

Because the number of dangling bonds is half of the number of three-fold-coordinated molecules, it follows that NF = NL

while for 4

For ESPP, the F- and L-type molecules may appear in one of the three edges consisting of hydrogen-bonded loops of four molecules, while in INTs, these molecules appear in the pentagons/squares, terminating the INT sprecies (see also Figure 1). Then, the following holds:

bLC + bCL = NL

i



(17)

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We wish to 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, CoE (2006-2011), and Lastu programs.

(18)



(19)

(20)

(21)

Another relation that holds in all studied clusters is (22)

and as a result of eqs 22−17 bFL + bLF = 2(NF − bFF)

REFERENCES

(1) Lin, S.-S. Detection of Large Water Clusters by a Low rf Quadrupole Mass Filter. Rev. Sci. Instrum. 1973, 44, 516−517. (2) Zwier, T. S. The Structure of Protonated Water Clusters. Science 2004, 304, 1119−1120. (3) Shi, Z.; Ford, J. V.; Wei, S.; Castleman, A. W., Jr. Water Clusters: Contributions of Binding Energy and Entropy to Stability. J. Chem. Phys. 1993, 99, 8009−8015. (4) Parkkinen, P.; Riikonen, S.; Halonen, L. Global Minima of Protonated Water Clusters (H2O)20H+ Revisited. J. Phys. Chem. A 2012, 116, 10826−10835. (5) Kuo, J.-L.; Klein, M. L. Structure of Protonated Water Clusters: Low-Energy Structures and Finite Temperature Behavior. J. Chem. Phys. 2005, 122, 024516. (6) Wu, C.-C.; Lin, C.-K.; Chang, H.-C.; Jiang, J.-C.; Kuo, J.-L.; Klein, M. L. Protonated Clathrate Cages Enclosing Neutral Water Molecules: H+(H2O)21 and H+(H2O)28. J. Chem. Phys. 2005, 122, 074315. (7) Laasonen, K.; Klein, M. L. Structural Study of (H2O)20 and (H2O)21H+ Using Density Functional Methods. J. Phys. Chem. 1994, 98, 10079−10083. (8) Hodges, M. P.; Wales, D. J. Global Minima of Protonated Water Clusters. Chem. Phys. Lett. 2000, 324, 279−288. (9) Arshad, K. 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. (10) 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. (11) James, T.; Wales, D. J. Protonated Water Clusters Described by an Empirical Valence Bond Potential. J. Chem. Phys. 2005, 122, 134306. (12) 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. (13) 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.

From eqs 18−20, it follows that

bLL = bFF

AUTHOR INFORMATION

Notes

bCF + bCL = 2NH2O − 4NF − bCC = 2NH2O − 4NF

bCL = bFC

(25)

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

Because the number of bonds that the C-type molecules donate is 2NC = 2(NH2O − NF − NL) = 2NH2O − 4NF, the following holds:

bLC = bCF

i

Corresponding Author

bCC = Nb − (bFF + bLL + bFL + bLF) − (bFC + bCF)

− (2NH2O − 5NF) = NF

∑ pi = 5

Both equations reduce the number of necessary fitting parameters.

A result of eqs 15−18 and the fact that none of the clusters contains anything but three- and four-fold-coordinated molecules gives

− (bLC + bCL) = 2NH2O − 5NF

4

∑ si = 5(m − 1)

In the prismatic and cubic morphologies considered in this article, the C-type molecules (i.e., four-fold-coordinated) reside near the center of mass of the cluster. Connecting two C-type molecules with one another via hydrogen bonds can be made using bonds between the C-type molecules only. Also, in these clusters, every three-fold-coordinated molecule is a neighbor of one C-type molecule. Then bFC + bCF = NF

clusters

s1 = −0.5(bFC + s3 + s2) + 2.5(m − 1)

(16)

bFF + bLL + bFL + bLF = NL + NF = 2 × NF

INTm5

(23)

The aTh chains appear naturally in the hydrogen-bonded loops that contain the F- and L-type molecules. In INTm4 structures, the numbers of F−F−F and L−L−L chains are equal (aFFF = aLLL) because if the four-molecule loops at one end of the cluster contain three F molecules, there must be a aTh chain. In this case, the other end has to contain an L−L−L chain because NF = NL. This does not hold for INTm5 ’s and ESPP as the rings at the different ends of a cluster can have different orders. Besides the dependencies mentioned above, some non-trivial dependencies with novel connectivity parameters were found during the fitting procedure, when unnatural bond angles were excluded (this imposes constraints on the hydrogen-bonded network). For INTm4 clusters, the following relations hold 9996

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

(14) 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. (15) Garczarek, F.; Gerwert, K. Functional Waters in Intraprotein Proton Transfer Monitored by FTIR Difference Spectroscopy. Nature 2006, 439, 109−112. (16) Niedner-Schatteburg, G.; Bondybey, V. E. FT-ICR Studies of Solvation Effects in Ionic Water Cluster Reactions. Chem. Rev. 2000, 100, 4059−4086. (17) 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. (18) Vaida, V. Perspective: Water Cluster Mediated Atmospheric Chemistry. J. Chem. Phys. 2011, 135, 020901. (19) 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. (20) 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. (21) Kazimirski, J. K.; Buch, V. Search for Low Energy Structures of Water Clusters (H2O)n, n = 20−22, 48, 123, and 293. J. Phys. Chem. A 2003, 107, 9762−9775. (22) Millot, C.; Soetens, J.-C.; Martins Costa, M. T. C.; Hodges, M. P.; Stone, A. J. Revised Anisotropic Site Potentials for the Water Dimer and Calculated Properties. J. Phys. Chem. A 1998, 102, 754− 770. (23) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91, 6269−6271. (24) Kabrede, H.; Hentschke, R. Global Minima of Water Clusters (H2O)N, N ≤ 25, Described by Three Empirical Potentials. J. Phys. Chem. B 2003, 107, 3914−3920. (25) Lenz, A.; Ojamäe, L. On the Stability of Dense Versus CageShaped Water Clusters: Quantumchemical Investigations of ZeroPoint Energies, Free Energies, Basis-Set Effects and IR Spectra of (H2O)12 and (H2O)20. Chem. Phys. Lett. 2006, 418, 361−367. (26) Fanourgakis, G. S.; Aprà, E.; Xantheas, S. S. High-Level Ab Initio Calculations for the Four Low-Lying 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. (27) Lenz, A.; Ojamäe, L. A Theoretical Study ofWater Equilibria: The Cluster Distribution Versus Temperature and Pressure for (H2O)n, n = 1−60, and Ice. J. Chem. Phys. 2009, 131, 134302. (28) Tokmachev, A. M.; Dronskowski, R. Hydrogen-Bond Networks in Finite Ice Nanotubes. J. Comput. Chem. 2011, 32, 99−105. (29) Nigra, P.; Kais, S. Pivot Method for Global Optimization: A Study of Water Clusters (H2O)N with 2 ≤ N ≤ 33. Chem. Phys. Lett. 1999, 305, 433−438. (30) Anick, D. J. Topology-Energy Relationships and Lowest Energy Configurations for Pentagonal Dodecahedral (H2O)20X Clusters, X = Empty, H2O, NH3, H3O+: The Importance of O-Topology. J. Chem. Phys. 2010, 132, 164311. (31) Sediki, A.; Lebsir, F.; Martiny, L.; Dauchez, M.; Krallafa, A. Ab Initio Investigation of the Topology and Properties of ThreeDimensional Clusters of Water (H2O)n. Food Chem. 2008, 106, 1476−1484 ; 4th International Workshop on Water in Foods. (32) Shi, Q.; Belair, S. D.; Francisco, J. S.; Kais, S. On the Interactions between Atmospheric Radicals and Cloud Droplets: A Molecular Picture of the Interface. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 9686− 9690. (33) 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.

(34) 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. (35) Anick, D. J. Zero Point Energy of Polyhedral Water Clusters. J. Phys. Chem. A 2005, 109, 5596−5601. (36) Radhakrishnan, T. P.; Herndon, W. C. Graph Theoretical Analysis of Water Clusters. J. Phys. Chem. 1991, 95, 10609−10617. (37) 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. (38) Lenz, A.; Ojamäe, L. Theoretical IR Spectra for Water Clusters (H2O)n (n = 6−22, 28, 30) and Identification of Spectral Contributions from Different H-Bond Conformations in Gaseous and Liquid Water. J. Phys. Chem. A 2006, 110, 13388−13393. (39) 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. (40) Anick, D. J. Polyhedral Water Clusters, I: Formal Consequences of the Ice Rules. J. Mol. Struct.: THEOCHEM 2002, 587, 87−96. (41) Anick, D. J. Polyhedral Water Clusters, II: Correlations of Connectivity Parameters with Electronic Energy and Hydrogen Bond Lengths. J. Mol. Struct.: THEOCHEM 2002, 587, 97−110. (42) Bahn, S. R.; Jacobsen, K. W. An Object-Oriented Scripting Interface to a Legacy Electronic Structure Code. Comput. Sci. Eng. 2002, 4, 56−66. (43) 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. (44) 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. (45) 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. (46) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (47) 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. (48) Tkatchenko, A.; Scheffler, M. Accurate Molecular van der Waals Interactions from Ground- State Electron Density and Free-Atom Reference Data. Phys. Rev. Lett. 2009, 102, 073005. (49) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (50) Cramer, C. J. Essentials of Computational Chemistry: Theories and Models; John Wiley: New York, 2004. (51) Jiang, J.-C.; Wang, Y.-S.; Chang, H.-C.; Lin, S. H.; Lee, Y. T.; Niedner-Schatteburg, G.; Chang, H.-C. Infrared Spectra of H+(H2O)5−8 Clusters: Evidence for Symmetric Proton Hydration. J. Am. Chem. Soc. 2000, 122, 1398−1410. (52) Miyazaki, M.; Fujii, A.; Ebata, T.; Mikami, N. Infrared Spectroscopic Evidence for Protonated Water Clusters Forming Nanoscale Cages. Science 2004, 304, 1134−1137. (53) 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. (54) Tokmachev, A. M.; Tchougréeff, A. L.; Dronskowski, R. Hydrogen-Bond Networks in Water Clusters (H2O)20: An Exhaustive Quantum-Chemical Analysis. ChemPhysChem 2010, 11, 384−388. (55) Kuo, J.; Ciobanu, C.; Ojamäe, L.; Shavitt, I.; Singer, S. Short HBonds and Spontaneous Selfdissociation in (H2O)20: Effects of HBond Topology. J. Chem. Phys. 2003, 118, 3583−3588. (56) Jurecka, P.; Sponer, J.; Cerny, J.; Hobza, P. Benchmark Database of Accurate (MP2 and CCSD(T) Complete Basis Set Limit) 9997

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998

The Journal of Physical Chemistry A

Article

Interaction Energies of Small Model Complexes, DNA Base Pairs, and Amino Acid Pairs. Phys. Chem. Chem. Phys. 2006, 8, 1985−1993. (57) Cybulski, S. M.; Lytle, M. L. The Origin of Deficiency of the Supermolecule Second-Order Møller−Plesset Approach for Evaluating Interaction Energies. J. Chem. Phys. 2007, 127, 141102. (58) Kelkkanen, A. K.; Lundqvist, B. I.; Nørskov, J. K. Density Functional for van der Waals Forces Accounts for Hydrogen Bond in Benchmark Set of Water Hexamers. J. Chem. Phys. 2009, 131, 046102. (59) Bates, D. M.; Tschumper, G. S. CCSD(T) Complete Basis Set Limit Relative Energies for Low-Lying Water Hexamer Structures. J. Phys. Chem. A 2009, 113, 3555−3559. (60) Wang, Y.; Babin, V.; Bowman, J. M.; Paesani, F. The Water Hexamer: Cage, Prism, or Both. Full Dimensional Quantum Simulations Say Both. J. Am. Chem. Soc. 2012, 134, 11116−11119. (61) Lenz, A.; Ojamäe, L. A Theoretical Study of Water Clusters: The Relation between Hydrogen Bond Topology and Interaction Energy from Quantum-Chemical Computations for Clusters with up to 22 Molecules. Phys. Chem. Chem. Phys. 2005, 7, 1905−1911. (62) 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.

9998

dx.doi.org/10.1021/jp4003092 | J. Phys. Chem. A 2013, 117, 9985−9998