- Email: [email protected]
Low temperature heat capacity of water clusters
Chemical Physics Letters 610–611 (2014) 369–374
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Low temperature heat capacity of water clusters Hongshan Chen a,∗ , Klavs Hansen b a b
College of Physics & Electronic Engineering, Northwest Normal University, Lanzhou 730070, China Department of Physics, University of Gothenburg, SE-41296 Gothenburg, Sweden
a r t i c l e
i n f o
Article history: Received 22 February 2014 In final form 28 July 2014 Available online 2 August 2014
a b s t r a c t Geometry optimization and vibrational frequency calculation are carried out at the MP2/6-31G(d,p) level for 35 low-energy isomers of (H2 O)n clusters in the size range n = 6–21. The heat capacities of the clusters are calculated using quantum statistical theories based on the harmonic approximation. The specific heat capacity increases with the cluster size but the difference diminishes gradually with increasing size. The heat capacities divided by the number of intermolecular vibrational modes are very close for all the clusters. The overall picture of the heat capacity of the clusters is bulk-like and it agrees well with the experimental results of size-selected clusters. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Water clusters are important objects of current experimental and theoretical studies because of their importance in understanding atmospheric chemistry, solution chemistry, and a large number of biochemical processes [1,2]. Substantial progress has been made in the study of the structure of water clusters. From a theoretical point of view, unbiased global minima searches based on different model potentials have been studied with the aim of understanding the characteristics of hydrogen bonds [3–5]. A number of ab initio calculations have also been carried out to investigate the strength of the hydrogen bonds and their cooperativity [6–13]. Experimentally, laser infrared and far-infrared vibration–rotation–tunneling spectroscopic studies demonstrated conclusively cyclic structures for n = 3–5 [14,15], and different isomers were observed for n = 6 [15–17]. Experiments involving pure water clusters as well as hydrated molecules were also carried out for n = 7–10 [18–21]. There now exists a general consensus about the global minima up to n = 8. Ab initio calculations also agree on the structure of the lowest energy isomers for n = 10, 12, 15, 16, and 20. The low-energy structures have been well established up to the size around 20 even though the global minima based on different models differ at some sizes. The structures and energetics of water clusters have been recently reviewed [22,23]. Relative to the extensive studies on the structure, investigations on the thermal behavior were limited, and we know little about the property of water clusters at finite temperature.
∗ Corresponding author. E-mail addresses: [email protected], [email protected] (H. Chen). http://dx.doi.org/10.1016/j.cplett.2014.07.063 0009-2614/© 2014 Elsevier B.V. All rights reserved.
Molecular dynamics and Monte Carlo simulations have been widely used to study the thermal properties of water clusters [24–34]. However, the reliability of the simulations depends on the model potentials used and the results can be different markedly. Moreover, the conventional MD and MC simulations were performed at the level of classical dynamics to treat the motion of atoms or rigid molecules. As the discrete vibration energy levels are allowed to change continuously, the calculations strongly overestimate the energy of the system and the heat capacities obtained from classical MD or MC simulations are expected to be much larger than the true value. We are aware of three literature studies that were carried out to simulate water clusters and treated atomic motion quantum mechanically [35–37], but only one was aimed to calculate the heat capacity. The variational gaussian wavepacket combined with a model potential SPC was applied to calculate the heat capacities of (H2 O)8 and (H2 O)10 [37]. The results show that the classical heat capacity is about one order larger than the quantum value. Simulations dealing with both the motion of the mass center and interatomic interactions in quantum mechanics are beyond the computation ability at present. Experimentally, the heat capacity of water cluster anions consisting of 48 and 118 molecules were measured using photofragmentation technique [38], and the protonated water clusters with 60–79 molecules were studied by multi-collision excitation of the accelerated clusters with helium [39]. While the phase transition temperatures are higher for the protonated clusters, the measured caloric curves of the anions and the protonated clusters coincide at the temperatures before the phase transition occurs and the overall picture is bulk-like. In this Letter, we turn to the quantum statistical theory to calculate the heat capacity of water clusters. Geometry optimization and vibrational frequencies were computed by the MP2/6-31G(d,p)
370
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
method for the low energy isomers of water clusters from size 6 to 21. The heat capacities of the clusters were obtained using the quantum statistical theories. The Letter was organized as follows. The computation methods are described in Section 2 and the main results reported in Section 3. After the discussion Section 4, the article ends with a brief summary. 2. Calculation methods In the harmonic approximation, the normal modes are independent and the vibrational motion is the superposition of the normal modes. At normal temperatures (say about 300 K), the modes with high frequencies remain frozen at the vibrational ground levels while the low frequency modes are partially populated at the excited states with probability exp(−h/kB T ). When temperature changes, the heat capacity contributed from the vibrational degrees of freedom will change dramatically. For a system having N normal modes of frequencies i (i = 1 − N), the energy states are E=
N
ni +
i=1
1 2
hi
(1)
For a canonical ensemble, the partition function is defined as [40] Z=
∞
e
−ˇ
N i=1
ni + 1 hi 2
=
ni =0
N e−ˇhi /2 i=1
(2)
1 − e−ˇhi
Here ˇ = 1/kB T . At temperature T, the average energy of the system is given by
hi ∂ ln Z = U0 + hi ⁄k T ∂ˇ B −1 i=1 e N
U=−
∂U ∂T
= kB V
N 2 h
e
i
i=1
kB T
(e
hi ⁄k T B
hi ⁄k T B − 1)
Clusters
Eb (eV) (without ZPE)
Eb (eV)
Eb /n
w6(prism) w6(cage) w6(ring) w7 w8 w9 w10 w10(3p) w11 w11(HF) w11(5p) w12 w12(3p) w12(5p) w13 w13(5p) w14 w14(HF) w14(5p) w15 w15(5p) w16 w16(5p) w17 w17(5p) w18 w18(5p) w19 w19(HF) w20 w20(4p) w20(5p) w21 w21(5p)
−2.89 −2.83 −2.57 −3.51 −4.43 −4.92 −5.60 −5.30 −6.18 −6.16 −6.05 −7.18 −6.82 −6.79 −7.68 −7.40 −8.44 −8.13 −8.03 −8.94 −8.75 −9.93 −9.26 −10.46 −10.22 −11.20 −10.77 −11.80 −11.71 −12.74 −12.60 −12.38 −13.24 −13.09
−2.22 −2.17 −1.98 −2.72 −3.47 −3.86 −4.40 −4.21 −4.86 −4.83 −4.76 −5.66 −5.37 −5.36 −6.06 −5.84 −6.68 −6.42 −6.35 −7.23 −7.11 −8.04 −7.56 −8.48 −8.31 −9.07 −8.77 −9.62 −9.50 −10.33 −10.23 −10.11 −10.77 −10.70
−0.37 −0.36 −0.33 −0.39 −0.43 −0.43 −0.44 −0.42 −0.44 −0.44 −0.43 −0.47 −0.45 −0.45 −0.47 −0.45 −0.48 −0.46 −0.45 −0.48 −0.47 −0.50 −0.47 −0.50 −0.49 −0.50 −0.49 −0.51 −0.50 −0.52 −0.51 −0.51 −0.51 −0.51
(3)
where U0 is the zero point energy. The heat capacity of the system can be derived as CV =
Table 1 Binding energy Eb and averaged over the number of molecules for (H2 O)n (n = 6–21).
2
(4)
The heat capacity is the sum of the individual oscillator heat capacities. Møller–Plesset 2nd order perturbation theory (MP2) was adopted to calculate the frequencies in our computation. The basis set used is the double split basis including polarization functions on both oxygen and hydrogen (6-31G(d,p)). The calculation was carried out using the gaussian 03 program [41]. Comparison with the experimental IR spectra of water clusters shows that the frequencies calculated at MP2/6-31G(d,p) level are quite satisfying [9,16,18,20]. The comparison of the calculated frequencies of the water octamer with the experimental IR spectra is given in supplementary material Figure S1. Experimental spectra are only available for the OH stretching band and the calculated frequencies are 2.0–5.6% higher. We also computed the vibrational frequencies of the water tetramer using a very large basis set 6-311++g(3df,3pd). The results are given in the supplementary material Table S1, and it shows the frequencies calculated at the two basis sets agree well. The water clusters considered in the present study are the sizes n = 6–21. The isomers of the clusters include the global minima structures based on empirical TIPnP (n = 3, 4, 5) potentials [3,5] and HF/6-31G(d,p) calculations [6]. (H2 O)6 has a few low energy isomers with competing energies. For n = 7–9, optimization based on different empirical potentials and ab initial calculations lead to same global minima structure. In the size range n = 10–21, the lowest energy structure based on different models differs amongst a few isomers. Totally we computed the frequencies of 35 structures.
The initial coordinate parameters of TIP4P/5P and HF isomers were taken from the Cambridge Cluster Database posted by Wales and his colleagues [42]. The TIP3P structures are based on our genetic algorithm optimization and they agree with those reported by Kabrede and Hentschke [4]. All the structures are illustrated in the supplementary material Figure S2. The isomers of water clusters in this Letter are denoted as following: Wn is the lowest energy structure of (H2 O)n based on the present calculation and Wn(kp) (k = 3, 4, 5) or Wn(HF) refers to the global minima isomers based on model potentials TIPkP (k = 3, 4, 5) or HF/6-31G(d,p) calculations, respectively. 3. Results 3.1. Structure and stability of (H2 O)n (n = 6–21) In order to study the correlation of the heat capacities with the stability of the clusters, the binding energies of the clusters are presented first but briefly. The binding energy is defined as the energy difference of the cluster at the equilibrium structure and the H2 O molecules building up it: Eb = E0 ((H2 O)n ) − n × E0 (H2 O)
(5)
The binding energies corrected and uncorrected with zero point energy (ZPE) are presented in Table 1. Ab initio calculations have been previously carried out for water hexamer [6–8,12], for the sizes n = 7–10 [6,10,12], n = 11–13[9,10,12,13], n = 14–16 [10,12], and n = 17–21 [10–12]. Even the binding energies calculated at different levels may be different slightly, the order of the stability of the isomers is in good agreement. Water hexamer is the smallest size that has 3D structures as global minimum isomer. Calculations at high levels show that the prism and cage forms are nearly isoenergetic [8]. The global minimum structure for n = 7, 8, 9 have
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
371
been well confirmed [9,10,15,19–21]. The lowest energy isomer of (H2 O)10 obtained by most model potentials is fused pentamers and the isomer based on TIP3P is a butterfly structure. The IR spectra seem to support the butterfly structure [20], but ab initio calculations show that the pentameric rings is lower in energy [10,12]. Calculations at MP4 level with the 6-311++G(2d,2p) basis set were employed to examine the structure of (H2 O)11 and (H2 O)13 [13]. The results suggest that there exist a few near-isoenergetic isomers for both clusters, but the lowest energy structures are in accordance with our calculation. Based on unbiased global minima optimization using empirical potentials and on ab initio calculations, the most stable structures are linear stacked cubes for n = 12, 16, 20 and stacked pentamers for n = 15 [6,10–12]. The stackedcube structures (n = 8, 12, 16, 20) are particularly more stable than their neighbors. n = 17–21 is the size range where the transition from all-surface cage to molecule-centered structure occurs. Based on the TTM2-F potential, the global minimum structure of (H2 O)17 has one interior molecule and it has been confirmed by MP2 calculation, but this structure is only 0.02 eV more stable than the lowest energy structure considered here. There are two other structures with one interior water molecule in this study. One is the lowest energy structure of (H2 O)19 and the other is the isomer of (H2 O)21 obtained by TIP5P. The binding energies of the all-surface cage and the molecule-centered isomer are very close in both cases. 3.2. Low temperature heat capacity of water clusters To compare the heat capacities of different clusters, the values are divided by the number of the water molecules in each cluster. The heat capacities per water molecule are reported in Figure 1. For reason of clarity, they are illustrated separately for the smaller and the larger clusters. The results of (H2 O)13 and (H2 O)14 are presented in both panels. It is obvious from the figure that heat capacities per molecule depend on both the size and the stability of clusters. Cv /n increases with the cluster size, but the difference diminishes gradually with increasing cluster size and the values saturate around size n = 13. The heat capacity is also correlated with the stability of the clusters. The heat capacity of the most stable isomer is smaller than that of other isomers at all the sizes. While the value of Cv /n of w16 (one of the most stable structure) is smaller than the values of the smaller clusters w13(5p), w14(HF),w14(5p), and w15(5p), Cv /n of the less stable structures w10(3p) and w11(5p) are greater than the values of the larger clusters w12, w13 and w14. The cyclic hexamer is the only ring structure studied here. Its heat capacity (per molecule) is much larger than that of other clusters at low temperature. We calculated the mean value of Cv /n and the root mean square deviation over the clusters with n = 13–21. The rms deviation is within 5% at temperatures above 50 K. To compare the heat capacities of different clusters, a more rational choice is to average the values over the normal vibration modes. The intramolecular vibrations do not contribute to the heat capacity (calculation using Eq. (4) shows that the contribution from intramolecular degrees of freedom is only 0.02% even at 200 K) and the number of normal modes for intermolecular vibration is Nf = 6n − 6. The heat capacities divided by Nf were reported in Figure 2. Again they were illustrated separately for the smaller and the larger clusters. Cv /Nf in Figure 2 shows no pronounced size dependence except for w6 and w8, of which the heat capacities are smaller at low temperature. The heat capacities divided by the number of intermolecular vibrational modes correlate very weakly with the stability of the clusters. The values of the most stable structures (w8, w12, w14, w16, w20) are a little smaller while the values of the less stable structures (w11(3p) and wn(5p) for n = 11, 13, 14, 15, 16, 18, 20) are slightly larger. Smaller heat capacities for the more stable structures are anticipated. Larger stability means stronger hydrogen bonds holding the water molecules together
Figure 1. Heat capacity Cv /n of water clusters by quantum statistics based on the harmonic frequencies calculated by MP2/6-31G(d,p). The experimental heat capacity of bulk ice is from Refs. [44,45].
and hence higher frequencies of the intermolecular vibrations, giving smaller heat capacities of the clusters. In fact, the results in Figure 2 show that the values of Cv /Nf are very close for all the clusters. The mean value of Cv /Nf and its root mean square deviation are calculated for all the clusters and the results are reported in Table 2. When the temperature increases, more vibrational modes are activated gradually and the heat capacity increases, but the rms deviation, even the absolute value, decreases. Table 2 Heat capacities of (H2 O)n averaged over the number of intermolecular vibrational modes. The root mean square deviations are presented in parenthesis. T (K)
Cv /Nf kB
T(K)
Cv /Nf kB
10 20 30 40 50 60 70 80 90 100
0.002 (0.001) 0.031 (0.007) 0.073 (0.010) 0.112 (0.012) 0.148 (0.013) 0.182 (0.013) 0.214 (0.012) 0.246 (0.011) 0.276 (0.010) 0.304 (0.009)
110 120 130 140 150 160 170 180 190 200
0.331 (0.008) 0.357 (0.007) 0.382 (0.007) 0.405 (0.007) 0.427 (0.007) 0.449 (0.007) 0.470 (0.007) 0.490 (0.007) 0.509 (0.007) 0.528 (0.007)
372
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
Figure 3. Comparison of the heat capacities of water cluster n = 48. The experimental heat capacity is from Ref. [38]. The predicted values are calculated using the average data in Table 2 and the contributions from the translational and rotational degrees of freedom.
Figure 2. Heat capacity Cv /Nf of water clusters by quantum statistics based on the harmonic frequencies by MP2/6-31G(d,p).
4. Discussion To our knowledge, there are only three experimental works that measure the heat capacity of water clusters. The heat capacities of water cluster anions consisting of 48 and 118 molecules were investigated using photofragmentation mass spectra [38]. The measured heat capacities are bulk-like at temperatures below about 100 K and an abrupt increase in the heat capacities is observed at 93 K and 118 K for size n = 48 and 118, respectively. The caloric curves of protonated water clusters H+ (H2 O)n with n = 60–79 were also studied recently [39]. The technique is based on multi-collision excitation of the accelerated clusters with helium. The measured caloric curves of different cluster sizes essentially retrace each other. The transition temperatures of the protonated clusters are higher (133 ± 6 K), but the curves resemble those of water cluster anions and fit the bulk ice curve. Based on the experimentally measured evaporation rate of water monomers from the clusters, we estimated the heat capacities of charged H+ (H2 O)n and X− (H2 O)n (X = O2 , CO3 , or NO3 ) in the size range n = 5–300 [43]. The temperature of the clusters is around 160 K, and the heat capacities of the clusters containing 10–19 molecules are about the values of bulk ice above 210 K. Using the average values in Table 2 (multiplying by Nf = 6n − 6 = 282) we can get the heat capacity of (H2 O)48
contributed from the vibrational degrees of freedom. The rotational constant B (Erot = J(J + 1)B) of (H2 O)20 is about 0.17 GHz and the corresponding temperature B/kB is only 0.007 K. (The rotational constant is inversely proportional to the moment of inertia and the corresponding temperature is lower for (H2 O)48 .) So we can take the heat capacity from the translational and rotational degrees of freedom as 3kB . Adding together the contributions from the vibrational, rotational, and translational motions we rebuilt the heat capacity of water cluster (H2 O)48 . The values are illustrated with the experimental results in Figure 3. The heat capacities of the protonated clusters compare well to those of the cluster anions, indicating that the charge has a very small influence on the heat capacity of the clusters. Figure 3 shows the calculated and the measured heat capacities of (H2 O)48 are close to identical at temperatures up to 70 K. The experimental heat capacity of bulk ice [44,45] is also included in Figure 1. The temperature dependence of the computed heat capacities of water clusters is very similar to that of bulk ice. The values of the clusters are a little smaller than that of the bulk, especially at higher temperatures. One possible reason for this deviation is the overestimation of the calculated frequencies. Due to the incompleteness of the basis set and electron correlation treatment that is inherent in quantum mechanical calculations, scaling factors (<1.0) are usually used to correct the calculated frequencies. The effect of the anharmonicity requires the calculated harmonic frequencies be further scaled down. The role of anharmonicity in small water clusters was recently examined using vibrational perturbation theory (VPT2) [46]. The effect of the anharmonicity on the heat capacity of ice was analyzed theoretically [47] and estimated using the density of vibrational states measured at 15 K [48]. Comparison between the measured heat capacities and the values calculated on the harmonic frequencies shows that anharmonic contribution appears at temperatures above 80 K and increases with elevating the temperature. We scaled our MP2/6-31G(d,p) harmonic frequencies by 0.9 and computed the heat capacities of the water clusters. The results show that the calculated heat capacities trace the experimentally measured curve of bulk ice up to 150 K but falls below at higher temperatures. As the predicted heat capacities and the experimental values of (H2 O)48 agree very well before the melting transition occurs, however, the difference between the heat capacities of the clusters and the bulk may have other reasons. For both bulk and clusters, the heat capacities are decided by the vibrational frequencies. We compute the densities of the
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
373
Figure 4. Densities of states of four phases of ice (left, from Ref. [48]) and four water clusters (right). The densities of the water clusters are normalized according to the number of the intermolecular vibration models. The values for ice lda, Ih and Ic (and for W15, W18 and W19) are shifted upwards by 0.04, 0.08, and 0.12, respectively.
vibrational states of the water clusters and compare them with the states of bulk ice in Figure 4. The densities of vibrational states of the cubic and hexagonal crystalline ice (Ic and Ih), and of the low-density and high-density amorphous ice (lda and hda) were determined by inelastic neutron scattering [48]. The four clusters are selected for their different structure features. W12 is stacked cubes and W15 is stacked pentamers. W18 is a double layer structure, and W19 has one internally solvated molecule. The frequencies of the intramolecular vibration in water clusters are very high (HOH bending around 1700 cm−1 and the OH bond stretching from 3500 to 4000 cm−1 ) and these vibrations are frozen at the ground states even at 1000 K. The densities of states (DOS) in Figure 4 include only the intermolecular vibration modes. The vibrational frequencies of the clusters are represented by gaussian distributions with a half width of 1 THz, and the densities are normalized according to the number of the intermolecular vibration modes Nf = 6n − 6. Similar to the densities of the bulk, the DOS of the clusters can be roughly divided into two sub-bands, corresponding to the translational and vibrational modes of the bulk. For the bulk, especially for the crystalline phases, there exists a clear separation of 5 THz between these two bands. For the clusters, the vibrational frequencies are obviously more extended and the librational modes shift to higher frequencies. The higher vibration frequencies reflect that water clusters form fewer but stronger hydrogen bonds since the molecules are arranged in an optimum manner in the lowenergy isomers. The stronger intermolecular bonding and higher vibrational frequencies might be the reason for the smaller heat capacities of the clusters. The medium strength of the hydrogen bond is responsible for the exceptional thermal properties of water. For the intermolecular vibrations in water clusters, the frequencies distribute from about 50 cm−1 to around 1300 cm−1 , and the corresponding activation temperatures (h/kB ) are from below 100 K to nearly 2000 K. So most of the vibrational excited states are only populated with very small probabilities. For bulk ice and water clusters, the heat capacities at 100 K are about 2kB per molecule, which means that on average only two of the nine vibrational degrees of freedom have contribution to the heat capacity and the remaining seven are frozen out. The classical dynamics allow the vibrational energies to change continuously and will therefore strongly overestimate the energy of the system. We know of three simulations carried out on water clusters applying quantum theory to describe the motion of the mass center (atoms or rigid molecules) [35–37]. The reweighted random series approach for stereographic path integral
Figure 5. Comparison of the heat capacities of water clusters based on classical dynamics, quantum variational gaussian wavepacket (VGW, Ref. [37]), and the quantum statistical methods.
Monte Carlo was used to simulate water clusters [36]. The melting temperature of (H2 O)n for n = 3–8 obtained in this study is close to 400 K, much higher than the results decided by other methods. The variational gaussian wavepacket (VGW) in combination with the replica-exchange Monte Carlo was applied to calculate the heat capacities of (H2 O)8 and (H2 O)10 [37]. A revised model potential of SPC was used to simulate rigid water molecules. The heat capacities computed by different methods are compared in Figure 5. It shows that the heat capacities obtained from VGW simulations are still larger than the values decided from the quantum statistics. The heat capacities determined by the classical dynamics, however, are about one order larger than the quantum values. The temperature dependence is also different for the heat capacities offered by the classical and quantum theories. 5. Conclusion Geometry optimization and frequency calculation were carried out at MP2/6-31G(d,p) level for the low energy isomers of water clusters in the size range n = 6–21. Heat capacities of water clusters were computed using quantum statistical theories upon the
374
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
harmonic approximation. The results show that the heat capacities per water molecule increase with the cluster size, but the values saturate around n = 13. The heat capacities correlate weakly with the stability of the clusters; the values are a little smaller for the isomers with greater stability. The heat capacities divided by the number of intermolecular vibrational modes Nf = 6n − 6 depend on the stability and the size very weakly and in fact they are very close for all the clusters. Mean values were calculated over all the structures for Cv /Nf and over the size range n = 13–21 for Cv /n. The relative rms deviations are only a few percent at temperatures above 50 K. The present calculation shows that the heat capacities of water clusters are bulk like, and the calculated values agree very well with the experimental results of cluster anions and protonated water clusters before the melting transition occurs. Acknowledgements Financial supports from the University of Gothenburg Nanoparticle Platform, from National Science Foundation of China (NFSC-11164024) and from Northwest Normal University (nwnukjcxgc-03-62). Martin Schmidt and Vladimir Mandelshtam are acknowledged for sending us their raw data. We thank Chalmers Center for Computational Science and Engineering (C3 SE) for generous amounts of computer time and Gansu Computation Center for technical supports. Chen gratefully acknowledges the support from China Scholarship Council. Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:10.1016/j.cplett.2014.07.063. References [1] [2] [3] [4] [5] [6] [7] [8]
B. Hartke, Angew. Chem. Int. Ed. 41 (2002) 1468. K. Muller-Dethlefs, P. Hobza, Chem. Rev. 100 (2000) 143. D.J. Wales, M.P. Hodges, Chem. Phys. Lett. 286 (1998) 65. H. Kabrede, R. Hentschke, J. Phys. Chem. B107 (2003) 3914. T. James, D.J. Wales, J. Hernández-Rojas, Chem. Phys. Lett. 415 (2005) 302. C. Lee, H. Chen, G. Fitzgerald, J. Chem. Phys. 102 (1995) 1266. J.K. Gregory, D.C. Clary, J. Phys. Chem. 100 (1996) 18014. J. Kim, K.S. Kim, J. Chem. Phys. 109 (1998) 5886.
[9] H.M. Lee, S.B. Suh, K.S. Kim, J. Chem. Phys. 114 (2001) 10749. [10] S. Maheshwary, N. Patel, N. Sathyamurthy, A.D. Kulkarni, S.R. Gadre, J. Phys. Chem. A 105 (2001) 10525. [11] G.S. Fanourgakis, E. Apra, S.S. Xantheas, J. Chem. Phys. 121 (2004) 2655. [12] A. Lenz, L. Ojamae, Phys. Chem. Chem. Phys. 7 (2005) 1905. [13] S. Bulusu, S. Yoo, E. Apra, S.S. Xantheas, X.C. Zeng, J. Phys. Chem. A 110 (2006) 11781. [14] N. Pugliano, R.J. Saykally, Science 257 (1992) 1937. [15] R.N. Pribble, T.S. Zwier, Science 265 (1994) 75. [16] K. Kim, K.D. Jordan, T.S. Zwier, J. Am. Chem. Soc. 116 (1994) 11568. [17] K. Liu, M.G. Brown, C. Carter, R.J. Saykally, J.K. Gregory, D.C. Clary, Nature 381 (1996) 501. [18] J. Brudermann, M. Melzer, U. Buck, J.K. Kazimirski, J. Sadlej, V. Bush, J. Chem. Phys. 110 (1999) 10649. [19] C.J. Gruenloh, J.R. Carney, C.A. Arrington, T.S. Zwier, S.Y. Fredericks, K.D. Jordan, Science 276 (1997) 1678. [20] U. Buck, I. Ettischer, M. Melzer, V. Buch, V. Sadlej, Phys. Rev. Lett. 80 (1998) 2578. [21] W.B. Blanton, et al., J. Am. Chem. Soc. 121 (1999) 3551. [22] S. Yoo, S.S. Xantheas, in: J. Leszczynski (Ed.), Handbook of Computational Chemistry, vol. 2, Springer, New York, 2012. [23] J.C. Howard, G.S. Tschumper, Wiley Interdisc. Rev.: Comput. Mol. Sci. 4 (2014) 169. [24] C.J. Tsai, K.D. Jordan, J. Chem. Phys. 99 (1993) 6957. [25] S.C. Farantos, S. Kapetanakis, A. Vegiri, J. Phys. Chem. 97 (1993) 12158. [26] M. Svanberg, J.B.C. Pettersson, J. Phys. Chem. A 102 (1998) 1865. [27] C.W. David, J. Chem. Phys. 73 (1980) 5395. [28] A. Vegiri, S.C. Farantos, J. Chem. Phys. 98 (1993) 4059. [29] P. Nigra, M.A. Carignano, S. Kaisa, J. Chem. Phys. 115 (2001) 2621. [30] A.N. Tharrington, K.D. Jordan, J. Phys. Chem. A 107 (2003) 7380. [31] D. Laria, J. Rodriguez, C. Dellago, D. Chandler, J. Phys. Chem. A 105 (2001) 2646. [32] M.A. Carignano, Chem. Phys. Lett. 361 (2002) 291. [33] A.V. Egorov, E.N. Brodskaya, A. Laaksonen, Synth. React. Inorg. Met.-Org. NanoMet. Chem. 38 (2008) 62. [34] T. Kaneko, T. Akimoto, K. Yasuoka, A. Mitsutake, X.C. Zeng, J. Chem. Theory Comput. 7 (2011) 3083. [35] S. Shin, W.S. Son, S. Jang, J. Mol. Struct.: THEOCHEM 673 (2004) 109. [36] S.F. Langley, E. Curotto, D.L. Freeman, J.D. Doll, J. Chem. Phys. 126 (2007) 084506. [37] P.A. Frantsuzov, V.A. Mandelshtam, J. Chem. Phys. 128 (2008) 94304. [38] C. Hock, et al., Phys. Rev. Lett. 103 (2009) 073401. [39] M. Schmidt, B.V. Issendorff, J. Chem. Phys. 136 (2012) 164307. [40] D.A. McQuarrie, Statistical Thermodynamics, Harper and Row, New York, 1973. [41] M.J. Frisch, et al., Gaussian 03, Revision D. 02, Gaussian, Inc., Wallingford, CT, 2004. [42] D.J. Wales et al., The Cambridge Cluster Database, http://www-wales.ch.cam. ac.uk/CCD.html [43] A.E.K. Sundén, et al., J. Chem. Phys. 130 (2009) 224308. [44] L.M. Shulman, Astron. Astrophys. 416 (2004) 187. [45] D.M. Murphy, T. Koop, Q. J. R. Meteorol. Soc. 131 (2005) 1539. [46] B. Temelso, G.C. Shields, J. Chem. Theory Comput. 7 (2011) 2804. [47] J. Leadbetter, Proc. R. Soc. London A 287 (1965) 403. [48] D.D. Klug, E. Whalley, E.C. Svensson, J.H. Root, V.F. Sears, Phys. Rev. B 44 (1991) 841.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Low temperature heat capacity of water clusters Hongshan Chen a,∗ , Klavs Hansen b a b
College of Physics & Electronic Engineering, Northwest Normal University, Lanzhou 730070, China Department of Physics, University of Gothenburg, SE-41296 Gothenburg, Sweden
a r t i c l e
i n f o
Article history: Received 22 February 2014 In final form 28 July 2014 Available online 2 August 2014
a b s t r a c t Geometry optimization and vibrational frequency calculation are carried out at the MP2/6-31G(d,p) level for 35 low-energy isomers of (H2 O)n clusters in the size range n = 6–21. The heat capacities of the clusters are calculated using quantum statistical theories based on the harmonic approximation. The specific heat capacity increases with the cluster size but the difference diminishes gradually with increasing size. The heat capacities divided by the number of intermolecular vibrational modes are very close for all the clusters. The overall picture of the heat capacity of the clusters is bulk-like and it agrees well with the experimental results of size-selected clusters. © 2014 Elsevier B.V. All rights reserved.
1. Introduction Water clusters are important objects of current experimental and theoretical studies because of their importance in understanding atmospheric chemistry, solution chemistry, and a large number of biochemical processes [1,2]. Substantial progress has been made in the study of the structure of water clusters. From a theoretical point of view, unbiased global minima searches based on different model potentials have been studied with the aim of understanding the characteristics of hydrogen bonds [3–5]. A number of ab initio calculations have also been carried out to investigate the strength of the hydrogen bonds and their cooperativity [6–13]. Experimentally, laser infrared and far-infrared vibration–rotation–tunneling spectroscopic studies demonstrated conclusively cyclic structures for n = 3–5 [14,15], and different isomers were observed for n = 6 [15–17]. Experiments involving pure water clusters as well as hydrated molecules were also carried out for n = 7–10 [18–21]. There now exists a general consensus about the global minima up to n = 8. Ab initio calculations also agree on the structure of the lowest energy isomers for n = 10, 12, 15, 16, and 20. The low-energy structures have been well established up to the size around 20 even though the global minima based on different models differ at some sizes. The structures and energetics of water clusters have been recently reviewed [22,23]. Relative to the extensive studies on the structure, investigations on the thermal behavior were limited, and we know little about the property of water clusters at finite temperature.
∗ Corresponding author. E-mail addresses: [email protected], [email protected] (H. Chen). http://dx.doi.org/10.1016/j.cplett.2014.07.063 0009-2614/© 2014 Elsevier B.V. All rights reserved.
Molecular dynamics and Monte Carlo simulations have been widely used to study the thermal properties of water clusters [24–34]. However, the reliability of the simulations depends on the model potentials used and the results can be different markedly. Moreover, the conventional MD and MC simulations were performed at the level of classical dynamics to treat the motion of atoms or rigid molecules. As the discrete vibration energy levels are allowed to change continuously, the calculations strongly overestimate the energy of the system and the heat capacities obtained from classical MD or MC simulations are expected to be much larger than the true value. We are aware of three literature studies that were carried out to simulate water clusters and treated atomic motion quantum mechanically [35–37], but only one was aimed to calculate the heat capacity. The variational gaussian wavepacket combined with a model potential SPC was applied to calculate the heat capacities of (H2 O)8 and (H2 O)10 [37]. The results show that the classical heat capacity is about one order larger than the quantum value. Simulations dealing with both the motion of the mass center and interatomic interactions in quantum mechanics are beyond the computation ability at present. Experimentally, the heat capacity of water cluster anions consisting of 48 and 118 molecules were measured using photofragmentation technique [38], and the protonated water clusters with 60–79 molecules were studied by multi-collision excitation of the accelerated clusters with helium [39]. While the phase transition temperatures are higher for the protonated clusters, the measured caloric curves of the anions and the protonated clusters coincide at the temperatures before the phase transition occurs and the overall picture is bulk-like. In this Letter, we turn to the quantum statistical theory to calculate the heat capacity of water clusters. Geometry optimization and vibrational frequencies were computed by the MP2/6-31G(d,p)
370
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
method for the low energy isomers of water clusters from size 6 to 21. The heat capacities of the clusters were obtained using the quantum statistical theories. The Letter was organized as follows. The computation methods are described in Section 2 and the main results reported in Section 3. After the discussion Section 4, the article ends with a brief summary. 2. Calculation methods In the harmonic approximation, the normal modes are independent and the vibrational motion is the superposition of the normal modes. At normal temperatures (say about 300 K), the modes with high frequencies remain frozen at the vibrational ground levels while the low frequency modes are partially populated at the excited states with probability exp(−h/kB T ). When temperature changes, the heat capacity contributed from the vibrational degrees of freedom will change dramatically. For a system having N normal modes of frequencies i (i = 1 − N), the energy states are E=
N
ni +
i=1
1 2
hi
(1)
For a canonical ensemble, the partition function is defined as [40] Z=
∞
e
−ˇ
N i=1
ni + 1 hi 2
=
ni =0
N e−ˇhi /2 i=1
(2)
1 − e−ˇhi
Here ˇ = 1/kB T . At temperature T, the average energy of the system is given by
hi ∂ ln Z = U0 + hi ⁄k T ∂ˇ B −1 i=1 e N
U=−
∂U ∂T
= kB V
N 2 h
e
i
i=1
kB T
(e
hi ⁄k T B
hi ⁄k T B − 1)
Clusters
Eb (eV) (without ZPE)
Eb (eV)
Eb /n
w6(prism) w6(cage) w6(ring) w7 w8 w9 w10 w10(3p) w11 w11(HF) w11(5p) w12 w12(3p) w12(5p) w13 w13(5p) w14 w14(HF) w14(5p) w15 w15(5p) w16 w16(5p) w17 w17(5p) w18 w18(5p) w19 w19(HF) w20 w20(4p) w20(5p) w21 w21(5p)
−2.89 −2.83 −2.57 −3.51 −4.43 −4.92 −5.60 −5.30 −6.18 −6.16 −6.05 −7.18 −6.82 −6.79 −7.68 −7.40 −8.44 −8.13 −8.03 −8.94 −8.75 −9.93 −9.26 −10.46 −10.22 −11.20 −10.77 −11.80 −11.71 −12.74 −12.60 −12.38 −13.24 −13.09
−2.22 −2.17 −1.98 −2.72 −3.47 −3.86 −4.40 −4.21 −4.86 −4.83 −4.76 −5.66 −5.37 −5.36 −6.06 −5.84 −6.68 −6.42 −6.35 −7.23 −7.11 −8.04 −7.56 −8.48 −8.31 −9.07 −8.77 −9.62 −9.50 −10.33 −10.23 −10.11 −10.77 −10.70
−0.37 −0.36 −0.33 −0.39 −0.43 −0.43 −0.44 −0.42 −0.44 −0.44 −0.43 −0.47 −0.45 −0.45 −0.47 −0.45 −0.48 −0.46 −0.45 −0.48 −0.47 −0.50 −0.47 −0.50 −0.49 −0.50 −0.49 −0.51 −0.50 −0.52 −0.51 −0.51 −0.51 −0.51
(3)
where U0 is the zero point energy. The heat capacity of the system can be derived as CV =
Table 1 Binding energy Eb and averaged over the number of molecules for (H2 O)n (n = 6–21).
2
(4)
The heat capacity is the sum of the individual oscillator heat capacities. Møller–Plesset 2nd order perturbation theory (MP2) was adopted to calculate the frequencies in our computation. The basis set used is the double split basis including polarization functions on both oxygen and hydrogen (6-31G(d,p)). The calculation was carried out using the gaussian 03 program [41]. Comparison with the experimental IR spectra of water clusters shows that the frequencies calculated at MP2/6-31G(d,p) level are quite satisfying [9,16,18,20]. The comparison of the calculated frequencies of the water octamer with the experimental IR spectra is given in supplementary material Figure S1. Experimental spectra are only available for the OH stretching band and the calculated frequencies are 2.0–5.6% higher. We also computed the vibrational frequencies of the water tetramer using a very large basis set 6-311++g(3df,3pd). The results are given in the supplementary material Table S1, and it shows the frequencies calculated at the two basis sets agree well. The water clusters considered in the present study are the sizes n = 6–21. The isomers of the clusters include the global minima structures based on empirical TIPnP (n = 3, 4, 5) potentials [3,5] and HF/6-31G(d,p) calculations [6]. (H2 O)6 has a few low energy isomers with competing energies. For n = 7–9, optimization based on different empirical potentials and ab initial calculations lead to same global minima structure. In the size range n = 10–21, the lowest energy structure based on different models differs amongst a few isomers. Totally we computed the frequencies of 35 structures.
The initial coordinate parameters of TIP4P/5P and HF isomers were taken from the Cambridge Cluster Database posted by Wales and his colleagues [42]. The TIP3P structures are based on our genetic algorithm optimization and they agree with those reported by Kabrede and Hentschke [4]. All the structures are illustrated in the supplementary material Figure S2. The isomers of water clusters in this Letter are denoted as following: Wn is the lowest energy structure of (H2 O)n based on the present calculation and Wn(kp) (k = 3, 4, 5) or Wn(HF) refers to the global minima isomers based on model potentials TIPkP (k = 3, 4, 5) or HF/6-31G(d,p) calculations, respectively. 3. Results 3.1. Structure and stability of (H2 O)n (n = 6–21) In order to study the correlation of the heat capacities with the stability of the clusters, the binding energies of the clusters are presented first but briefly. The binding energy is defined as the energy difference of the cluster at the equilibrium structure and the H2 O molecules building up it: Eb = E0 ((H2 O)n ) − n × E0 (H2 O)
(5)
The binding energies corrected and uncorrected with zero point energy (ZPE) are presented in Table 1. Ab initio calculations have been previously carried out for water hexamer [6–8,12], for the sizes n = 7–10 [6,10,12], n = 11–13[9,10,12,13], n = 14–16 [10,12], and n = 17–21 [10–12]. Even the binding energies calculated at different levels may be different slightly, the order of the stability of the isomers is in good agreement. Water hexamer is the smallest size that has 3D structures as global minimum isomer. Calculations at high levels show that the prism and cage forms are nearly isoenergetic [8]. The global minimum structure for n = 7, 8, 9 have
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
371
been well confirmed [9,10,15,19–21]. The lowest energy isomer of (H2 O)10 obtained by most model potentials is fused pentamers and the isomer based on TIP3P is a butterfly structure. The IR spectra seem to support the butterfly structure [20], but ab initio calculations show that the pentameric rings is lower in energy [10,12]. Calculations at MP4 level with the 6-311++G(2d,2p) basis set were employed to examine the structure of (H2 O)11 and (H2 O)13 [13]. The results suggest that there exist a few near-isoenergetic isomers for both clusters, but the lowest energy structures are in accordance with our calculation. Based on unbiased global minima optimization using empirical potentials and on ab initio calculations, the most stable structures are linear stacked cubes for n = 12, 16, 20 and stacked pentamers for n = 15 [6,10–12]. The stackedcube structures (n = 8, 12, 16, 20) are particularly more stable than their neighbors. n = 17–21 is the size range where the transition from all-surface cage to molecule-centered structure occurs. Based on the TTM2-F potential, the global minimum structure of (H2 O)17 has one interior molecule and it has been confirmed by MP2 calculation, but this structure is only 0.02 eV more stable than the lowest energy structure considered here. There are two other structures with one interior water molecule in this study. One is the lowest energy structure of (H2 O)19 and the other is the isomer of (H2 O)21 obtained by TIP5P. The binding energies of the all-surface cage and the molecule-centered isomer are very close in both cases. 3.2. Low temperature heat capacity of water clusters To compare the heat capacities of different clusters, the values are divided by the number of the water molecules in each cluster. The heat capacities per water molecule are reported in Figure 1. For reason of clarity, they are illustrated separately for the smaller and the larger clusters. The results of (H2 O)13 and (H2 O)14 are presented in both panels. It is obvious from the figure that heat capacities per molecule depend on both the size and the stability of clusters. Cv /n increases with the cluster size, but the difference diminishes gradually with increasing cluster size and the values saturate around size n = 13. The heat capacity is also correlated with the stability of the clusters. The heat capacity of the most stable isomer is smaller than that of other isomers at all the sizes. While the value of Cv /n of w16 (one of the most stable structure) is smaller than the values of the smaller clusters w13(5p), w14(HF),w14(5p), and w15(5p), Cv /n of the less stable structures w10(3p) and w11(5p) are greater than the values of the larger clusters w12, w13 and w14. The cyclic hexamer is the only ring structure studied here. Its heat capacity (per molecule) is much larger than that of other clusters at low temperature. We calculated the mean value of Cv /n and the root mean square deviation over the clusters with n = 13–21. The rms deviation is within 5% at temperatures above 50 K. To compare the heat capacities of different clusters, a more rational choice is to average the values over the normal vibration modes. The intramolecular vibrations do not contribute to the heat capacity (calculation using Eq. (4) shows that the contribution from intramolecular degrees of freedom is only 0.02% even at 200 K) and the number of normal modes for intermolecular vibration is Nf = 6n − 6. The heat capacities divided by Nf were reported in Figure 2. Again they were illustrated separately for the smaller and the larger clusters. Cv /Nf in Figure 2 shows no pronounced size dependence except for w6 and w8, of which the heat capacities are smaller at low temperature. The heat capacities divided by the number of intermolecular vibrational modes correlate very weakly with the stability of the clusters. The values of the most stable structures (w8, w12, w14, w16, w20) are a little smaller while the values of the less stable structures (w11(3p) and wn(5p) for n = 11, 13, 14, 15, 16, 18, 20) are slightly larger. Smaller heat capacities for the more stable structures are anticipated. Larger stability means stronger hydrogen bonds holding the water molecules together
Figure 1. Heat capacity Cv /n of water clusters by quantum statistics based on the harmonic frequencies calculated by MP2/6-31G(d,p). The experimental heat capacity of bulk ice is from Refs. [44,45].
and hence higher frequencies of the intermolecular vibrations, giving smaller heat capacities of the clusters. In fact, the results in Figure 2 show that the values of Cv /Nf are very close for all the clusters. The mean value of Cv /Nf and its root mean square deviation are calculated for all the clusters and the results are reported in Table 2. When the temperature increases, more vibrational modes are activated gradually and the heat capacity increases, but the rms deviation, even the absolute value, decreases. Table 2 Heat capacities of (H2 O)n averaged over the number of intermolecular vibrational modes. The root mean square deviations are presented in parenthesis. T (K)
Cv /Nf kB
T(K)
Cv /Nf kB
10 20 30 40 50 60 70 80 90 100
0.002 (0.001) 0.031 (0.007) 0.073 (0.010) 0.112 (0.012) 0.148 (0.013) 0.182 (0.013) 0.214 (0.012) 0.246 (0.011) 0.276 (0.010) 0.304 (0.009)
110 120 130 140 150 160 170 180 190 200
0.331 (0.008) 0.357 (0.007) 0.382 (0.007) 0.405 (0.007) 0.427 (0.007) 0.449 (0.007) 0.470 (0.007) 0.490 (0.007) 0.509 (0.007) 0.528 (0.007)
372
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
Figure 3. Comparison of the heat capacities of water cluster n = 48. The experimental heat capacity is from Ref. [38]. The predicted values are calculated using the average data in Table 2 and the contributions from the translational and rotational degrees of freedom.
Figure 2. Heat capacity Cv /Nf of water clusters by quantum statistics based on the harmonic frequencies by MP2/6-31G(d,p).
4. Discussion To our knowledge, there are only three experimental works that measure the heat capacity of water clusters. The heat capacities of water cluster anions consisting of 48 and 118 molecules were investigated using photofragmentation mass spectra [38]. The measured heat capacities are bulk-like at temperatures below about 100 K and an abrupt increase in the heat capacities is observed at 93 K and 118 K for size n = 48 and 118, respectively. The caloric curves of protonated water clusters H+ (H2 O)n with n = 60–79 were also studied recently [39]. The technique is based on multi-collision excitation of the accelerated clusters with helium. The measured caloric curves of different cluster sizes essentially retrace each other. The transition temperatures of the protonated clusters are higher (133 ± 6 K), but the curves resemble those of water cluster anions and fit the bulk ice curve. Based on the experimentally measured evaporation rate of water monomers from the clusters, we estimated the heat capacities of charged H+ (H2 O)n and X− (H2 O)n (X = O2 , CO3 , or NO3 ) in the size range n = 5–300 [43]. The temperature of the clusters is around 160 K, and the heat capacities of the clusters containing 10–19 molecules are about the values of bulk ice above 210 K. Using the average values in Table 2 (multiplying by Nf = 6n − 6 = 282) we can get the heat capacity of (H2 O)48
contributed from the vibrational degrees of freedom. The rotational constant B (Erot = J(J + 1)B) of (H2 O)20 is about 0.17 GHz and the corresponding temperature B/kB is only 0.007 K. (The rotational constant is inversely proportional to the moment of inertia and the corresponding temperature is lower for (H2 O)48 .) So we can take the heat capacity from the translational and rotational degrees of freedom as 3kB . Adding together the contributions from the vibrational, rotational, and translational motions we rebuilt the heat capacity of water cluster (H2 O)48 . The values are illustrated with the experimental results in Figure 3. The heat capacities of the protonated clusters compare well to those of the cluster anions, indicating that the charge has a very small influence on the heat capacity of the clusters. Figure 3 shows the calculated and the measured heat capacities of (H2 O)48 are close to identical at temperatures up to 70 K. The experimental heat capacity of bulk ice [44,45] is also included in Figure 1. The temperature dependence of the computed heat capacities of water clusters is very similar to that of bulk ice. The values of the clusters are a little smaller than that of the bulk, especially at higher temperatures. One possible reason for this deviation is the overestimation of the calculated frequencies. Due to the incompleteness of the basis set and electron correlation treatment that is inherent in quantum mechanical calculations, scaling factors (<1.0) are usually used to correct the calculated frequencies. The effect of the anharmonicity requires the calculated harmonic frequencies be further scaled down. The role of anharmonicity in small water clusters was recently examined using vibrational perturbation theory (VPT2) [46]. The effect of the anharmonicity on the heat capacity of ice was analyzed theoretically [47] and estimated using the density of vibrational states measured at 15 K [48]. Comparison between the measured heat capacities and the values calculated on the harmonic frequencies shows that anharmonic contribution appears at temperatures above 80 K and increases with elevating the temperature. We scaled our MP2/6-31G(d,p) harmonic frequencies by 0.9 and computed the heat capacities of the water clusters. The results show that the calculated heat capacities trace the experimentally measured curve of bulk ice up to 150 K but falls below at higher temperatures. As the predicted heat capacities and the experimental values of (H2 O)48 agree very well before the melting transition occurs, however, the difference between the heat capacities of the clusters and the bulk may have other reasons. For both bulk and clusters, the heat capacities are decided by the vibrational frequencies. We compute the densities of the
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
373
Figure 4. Densities of states of four phases of ice (left, from Ref. [48]) and four water clusters (right). The densities of the water clusters are normalized according to the number of the intermolecular vibration models. The values for ice lda, Ih and Ic (and for W15, W18 and W19) are shifted upwards by 0.04, 0.08, and 0.12, respectively.
vibrational states of the water clusters and compare them with the states of bulk ice in Figure 4. The densities of vibrational states of the cubic and hexagonal crystalline ice (Ic and Ih), and of the low-density and high-density amorphous ice (lda and hda) were determined by inelastic neutron scattering [48]. The four clusters are selected for their different structure features. W12 is stacked cubes and W15 is stacked pentamers. W18 is a double layer structure, and W19 has one internally solvated molecule. The frequencies of the intramolecular vibration in water clusters are very high (HOH bending around 1700 cm−1 and the OH bond stretching from 3500 to 4000 cm−1 ) and these vibrations are frozen at the ground states even at 1000 K. The densities of states (DOS) in Figure 4 include only the intermolecular vibration modes. The vibrational frequencies of the clusters are represented by gaussian distributions with a half width of 1 THz, and the densities are normalized according to the number of the intermolecular vibration modes Nf = 6n − 6. Similar to the densities of the bulk, the DOS of the clusters can be roughly divided into two sub-bands, corresponding to the translational and vibrational modes of the bulk. For the bulk, especially for the crystalline phases, there exists a clear separation of 5 THz between these two bands. For the clusters, the vibrational frequencies are obviously more extended and the librational modes shift to higher frequencies. The higher vibration frequencies reflect that water clusters form fewer but stronger hydrogen bonds since the molecules are arranged in an optimum manner in the lowenergy isomers. The stronger intermolecular bonding and higher vibrational frequencies might be the reason for the smaller heat capacities of the clusters. The medium strength of the hydrogen bond is responsible for the exceptional thermal properties of water. For the intermolecular vibrations in water clusters, the frequencies distribute from about 50 cm−1 to around 1300 cm−1 , and the corresponding activation temperatures (h/kB ) are from below 100 K to nearly 2000 K. So most of the vibrational excited states are only populated with very small probabilities. For bulk ice and water clusters, the heat capacities at 100 K are about 2kB per molecule, which means that on average only two of the nine vibrational degrees of freedom have contribution to the heat capacity and the remaining seven are frozen out. The classical dynamics allow the vibrational energies to change continuously and will therefore strongly overestimate the energy of the system. We know of three simulations carried out on water clusters applying quantum theory to describe the motion of the mass center (atoms or rigid molecules) [35–37]. The reweighted random series approach for stereographic path integral
Figure 5. Comparison of the heat capacities of water clusters based on classical dynamics, quantum variational gaussian wavepacket (VGW, Ref. [37]), and the quantum statistical methods.
Monte Carlo was used to simulate water clusters [36]. The melting temperature of (H2 O)n for n = 3–8 obtained in this study is close to 400 K, much higher than the results decided by other methods. The variational gaussian wavepacket (VGW) in combination with the replica-exchange Monte Carlo was applied to calculate the heat capacities of (H2 O)8 and (H2 O)10 [37]. A revised model potential of SPC was used to simulate rigid water molecules. The heat capacities computed by different methods are compared in Figure 5. It shows that the heat capacities obtained from VGW simulations are still larger than the values decided from the quantum statistics. The heat capacities determined by the classical dynamics, however, are about one order larger than the quantum values. The temperature dependence is also different for the heat capacities offered by the classical and quantum theories. 5. Conclusion Geometry optimization and frequency calculation were carried out at MP2/6-31G(d,p) level for the low energy isomers of water clusters in the size range n = 6–21. Heat capacities of water clusters were computed using quantum statistical theories upon the
374
H. Chen, K. Hansen / Chemical Physics Letters 610–611 (2014) 369–374
harmonic approximation. The results show that the heat capacities per water molecule increase with the cluster size, but the values saturate around n = 13. The heat capacities correlate weakly with the stability of the clusters; the values are a little smaller for the isomers with greater stability. The heat capacities divided by the number of intermolecular vibrational modes Nf = 6n − 6 depend on the stability and the size very weakly and in fact they are very close for all the clusters. Mean values were calculated over all the structures for Cv /Nf and over the size range n = 13–21 for Cv /n. The relative rms deviations are only a few percent at temperatures above 50 K. The present calculation shows that the heat capacities of water clusters are bulk like, and the calculated values agree very well with the experimental results of cluster anions and protonated water clusters before the melting transition occurs. Acknowledgements Financial supports from the University of Gothenburg Nanoparticle Platform, from National Science Foundation of China (NFSC-11164024) and from Northwest Normal University (nwnukjcxgc-03-62). Martin Schmidt and Vladimir Mandelshtam are acknowledged for sending us their raw data. We thank Chalmers Center for Computational Science and Engineering (C3 SE) for generous amounts of computer time and Gansu Computation Center for technical supports. Chen gratefully acknowledges the support from China Scholarship Council. Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:10.1016/j.cplett.2014.07.063. References [1] [2] [3] [4] [5] [6] [7] [8]
B. Hartke, Angew. Chem. Int. Ed. 41 (2002) 1468. K. Muller-Dethlefs, P. Hobza, Chem. Rev. 100 (2000) 143. D.J. Wales, M.P. Hodges, Chem. Phys. Lett. 286 (1998) 65. H. Kabrede, R. Hentschke, J. Phys. Chem. B107 (2003) 3914. T. James, D.J. Wales, J. Hernández-Rojas, Chem. Phys. Lett. 415 (2005) 302. C. Lee, H. Chen, G. Fitzgerald, J. Chem. Phys. 102 (1995) 1266. J.K. Gregory, D.C. Clary, J. Phys. Chem. 100 (1996) 18014. J. Kim, K.S. Kim, J. Chem. Phys. 109 (1998) 5886.
[9] H.M. Lee, S.B. Suh, K.S. Kim, J. Chem. Phys. 114 (2001) 10749. [10] S. Maheshwary, N. Patel, N. Sathyamurthy, A.D. Kulkarni, S.R. Gadre, J. Phys. Chem. A 105 (2001) 10525. [11] G.S. Fanourgakis, E. Apra, S.S. Xantheas, J. Chem. Phys. 121 (2004) 2655. [12] A. Lenz, L. Ojamae, Phys. Chem. Chem. Phys. 7 (2005) 1905. [13] S. Bulusu, S. Yoo, E. Apra, S.S. Xantheas, X.C. Zeng, J. Phys. Chem. A 110 (2006) 11781. [14] N. Pugliano, R.J. Saykally, Science 257 (1992) 1937. [15] R.N. Pribble, T.S. Zwier, Science 265 (1994) 75. [16] K. Kim, K.D. Jordan, T.S. Zwier, J. Am. Chem. Soc. 116 (1994) 11568. [17] K. Liu, M.G. Brown, C. Carter, R.J. Saykally, J.K. Gregory, D.C. Clary, Nature 381 (1996) 501. [18] J. Brudermann, M. Melzer, U. Buck, J.K. Kazimirski, J. Sadlej, V. Bush, J. Chem. Phys. 110 (1999) 10649. [19] C.J. Gruenloh, J.R. Carney, C.A. Arrington, T.S. Zwier, S.Y. Fredericks, K.D. Jordan, Science 276 (1997) 1678. [20] U. Buck, I. Ettischer, M. Melzer, V. Buch, V. Sadlej, Phys. Rev. Lett. 80 (1998) 2578. [21] W.B. Blanton, et al., J. Am. Chem. Soc. 121 (1999) 3551. [22] S. Yoo, S.S. Xantheas, in: J. Leszczynski (Ed.), Handbook of Computational Chemistry, vol. 2, Springer, New York, 2012. [23] J.C. Howard, G.S. Tschumper, Wiley Interdisc. Rev.: Comput. Mol. Sci. 4 (2014) 169. [24] C.J. Tsai, K.D. Jordan, J. Chem. Phys. 99 (1993) 6957. [25] S.C. Farantos, S. Kapetanakis, A. Vegiri, J. Phys. Chem. 97 (1993) 12158. [26] M. Svanberg, J.B.C. Pettersson, J. Phys. Chem. A 102 (1998) 1865. [27] C.W. David, J. Chem. Phys. 73 (1980) 5395. [28] A. Vegiri, S.C. Farantos, J. Chem. Phys. 98 (1993) 4059. [29] P. Nigra, M.A. Carignano, S. Kaisa, J. Chem. Phys. 115 (2001) 2621. [30] A.N. Tharrington, K.D. Jordan, J. Phys. Chem. A 107 (2003) 7380. [31] D. Laria, J. Rodriguez, C. Dellago, D. Chandler, J. Phys. Chem. A 105 (2001) 2646. [32] M.A. Carignano, Chem. Phys. Lett. 361 (2002) 291. [33] A.V. Egorov, E.N. Brodskaya, A. Laaksonen, Synth. React. Inorg. Met.-Org. NanoMet. Chem. 38 (2008) 62. [34] T. Kaneko, T. Akimoto, K. Yasuoka, A. Mitsutake, X.C. Zeng, J. Chem. Theory Comput. 7 (2011) 3083. [35] S. Shin, W.S. Son, S. Jang, J. Mol. Struct.: THEOCHEM 673 (2004) 109. [36] S.F. Langley, E. Curotto, D.L. Freeman, J.D. Doll, J. Chem. Phys. 126 (2007) 084506. [37] P.A. Frantsuzov, V.A. Mandelshtam, J. Chem. Phys. 128 (2008) 94304. [38] C. Hock, et al., Phys. Rev. Lett. 103 (2009) 073401. [39] M. Schmidt, B.V. Issendorff, J. Chem. Phys. 136 (2012) 164307. [40] D.A. McQuarrie, Statistical Thermodynamics, Harper and Row, New York, 1973. [41] M.J. Frisch, et al., Gaussian 03, Revision D. 02, Gaussian, Inc., Wallingford, CT, 2004. [42] D.J. Wales et al., The Cambridge Cluster Database, http://www-wales.ch.cam. ac.uk/CCD.html [43] A.E.K. Sundén, et al., J. Chem. Phys. 130 (2009) 224308. [44] L.M. Shulman, Astron. Astrophys. 416 (2004) 187. [45] D.M. Murphy, T. Koop, Q. J. R. Meteorol. Soc. 131 (2005) 1539. [46] B. Temelso, G.C. Shields, J. Chem. Theory Comput. 7 (2011) 2804. [47] J. Leadbetter, Proc. R. Soc. London A 287 (1965) 403. [48] D.D. Klug, E. Whalley, E.C. Svensson, J.H. Root, V.F. Sears, Phys. Rev. B 44 (1991) 841.