Low-Lying Structures and Stabilities of Large Water Clusters

Jul 29, 2010 - Appraisal of molecular tailoring approach for large clusters. Nityananda Sahu , Sachin D. Yeole , Shridhar R. Gadre. The Journal of Che...
0 downloads 11 Views 4MB Size
Download PDF
J. Phys. Chem. A 2010, 114, 9253–9261

9253

Low-Lying Structures and Stabilities of Large Water Clusters: Investigation Based on the Combination of the AMOEBA Potential and Generalized Energy-Based Fragmentation Approach Zhen Yang, Shugui Hua, Weijie Hua, and Shuhua Li* School of Chemistry and Chemical Engineering, Key Laboratory of Mesoscopic Chemistry of the Ministry of Education, Institute of Theoretical and Computational Chemistry, Nanjing UniVersity, Nanjing 210093, People’s Republic of China ReceiVed: April 28, 2010; ReVised Manuscript ReceiVed: July 15, 2010

Based on a large database of local minima obtained with the polarizable AMOEBA potential, the generalized energy-based fragmentation (GEBF) approach is applied to locate low-lying structures of water clusters (H2O)n in the range n ) 20-30, at the B3LYP and MP2 levels. Our results show that the relative stabilities of isomers predicted by the AMOEBA empirical potential differ noticeably from those predicted by GEBFB3LYP/6-311++G(d,p) and GEBF-MP2/6-311++G(3df,2p) calculations. From GEBF-B3LYP energies with zero-point vibrational energy corrections, one can see that for water clusters in the range n ) 20-30 the transition from one-centered to two-centered cage structure occurs at n ) 26. With increasing cluster size, the number of H-bonds per water molecule in the lowest-energy structures shows a gradually increasing trend, and the proportion of four-coordinated water molecules gradually increases, as expected for large water clusters. Based on GEBF-MP2/6-311++G(3df,2p) energies (instead of GEBF-B3LYP/6-311++G(d,p) energies), different lowest-energy structures can be found for six cluster sizes in the range n ) 20-30, suggesting the significance of the dispersion interaction in determining the relative energies of low-lying water clusters. 1. Introduction Understanding and predicting the properties of water clusters are of fundamental importance in many physical, chemical, and biological fields. An important purpose for such studies is to obtain a molecular-level description of the properties of the bulk water through the study of various water clusters, (H2O)n.1-3 Recently, a large number of experimental3-14 and theoretical studies15-43 have been conducted to investigate the structures and properties of water clusters. Since the number of local minimum structures for (H2O)n increases exponentially with the size of water clusters, many effective optimization approaches have been developed to search the global minima of water clusters with different sizes by using empirical interaction potentials.15-25 For example, Niesse and Mayne15 used a genetic algorithm with the TIP3P44 potential to explore the global minima of water clusters containing up to 13 molecules. Wales and co-workers18,19 employed a basin-hopping (BH) method to find the global minima for n e 21 with the TIP4P44 and TIP5P45 potentials. Kazachenko and Thakkar25 developed an improved minima-hopping algorithm to locate the global minima for n e 37. On the other hand, one can notice that different empirical potentials usually lead to different global minima, especially for relatively large water clusters (n > 10). For example, the global minima predicted by the TTM2-F potential agree qualitatively with those predicted by the TIP4P potential for small clusters n ) 2-11, but the results for larger clusters are not consistent with each other.21 The global minimum structures for n ) 11-21 predicted with TIP4P and TIP5P potentials are * To whom correspondence should be addressed. E-mail: shuhua@ nju.edu.cn.

quite different.19 Due to different accuracies of various empirical potentials, global minimum structures predicted from various model potentials are better considered as good candidates for true global minima. Thus, it is still necessary to reoptimize the structures of these low-energy candidates using high-level ab initio calculations. Presently, a number of ab initio calculations have been performed for the low-lying structures of water clusters.26-43 For example, Maheshwary et al.32 have extensively investigated structures and stabilities of water clusters with n ) 8-20 at the Hartree-Fock (HF) levels with up to the 6-311++G(2d,2p) basis set. Fanourgakis et al.33 reported more accurate energies for four low-lying families of (H2O)20 at the second-order Møller-Plesset perturbation (MP2) level of theory. Lenz and Ojamae35 performed geometry optimizations for a variety of water clusters with 8, 10, and 12 molecules, as well as for some larger clusters with up to 22 molecules, whose structures were constructed based on the information from small clusters. Up to now, the lowest-energy structures for n e 10 have been well established from various ab initio calculations, which are consistent with experimental studies.26-30,39-43 However, as mentioned above, ab initio calculations on larger water clusters (n g 20) have been limited to a relatively small number of low-energy structures.32-37 It should be emphasized that an exhaustive search for the global minimum structures of large water clusters is not yet practical for current ab initio calculations. Evidently, a feasible way of determining low-energy structures of water clusters is to first identify a large number of low-energy candidate structures with empirical potentials, and to further optimize the structures of these candidates using various ab initio methods. It is well-known that the computational cost of traditional algorithms for ab initio calculations increases steeply with the

10.1021/jp1038267  2010 American Chemical Society Published on Web 07/29/2010

9254

J. Phys. Chem. A, Vol. 114, No. 34, 2010

system size. To speed up ab initio calculations for large molecules, various linear-scaling algorithms have been developed.46-83 Among them, energy-based fragmentation approaches67-83 have been proven to be practical tools that can fairly accurately reproduce traditional ab initio results for large closed-shell systems. Especially the generalized energy-based fragmentation (GEBF) approach developed by our group71 has been demonstrated to be suitable for systems with many polar (or charged) groups such as water clusters.71-74 In the GEBF approach, the energy of a large molecule is directly evaluated as the linear combination of many small subsystems, which are embedded in the presence of background point charges. This approach can be applicable for various ab initio methods such as HF, MP2, density functional theory (DFT), and so on. Furthermore, this approach can be easily implemented with existing quantum chemistry packages, and highly efficient parallel computations can be achieved in the GEBF calculations. In the present work, our aim is to systematically investigate low-energy structures and their stabilities for water clusters with n ) 20-30 by using the combination of the empirical potential and GEBF approach. First, the potential energy surfaces (PESs) of water clusters are explored by using a modified BH method with the polarizable AMOEBA84 empirical potential. Then we carry out a series of GEBF-DFT calculations to determine lowenergy structures from a large number of local minimum structures obtained with the AMOEBA potential. The main purpose of this study is to obtain qualitative structural features for low-energy water clusters with different sizes (n ) 20-30), rather than to find true global minima for these clusters. This is mainly because an exhaustive sampling for the PESs of such large clusters with the AMOEBA potential would be computationally demanding, and our sampling in this study may be not sufficient. To the best of our knowledge, nevertheless, this study is the first attempt to systematically search low-energy minima of large water clusters (n ) 20-30) with a combination of empirical potential and ab initio calculations. The paper is arranged as follows. In section 2, we will first describe the modified BH optimization method and the AMOEBA potential, and then introduce the GEBF approach and a hybrid procedure for locating low-energy structures. In section 3, we will give 10 lowest-energy structures for each water cluster obtained with the AMOEBA potential and with ab initio calculations, respectively. Their qualitative structural features will be discussed and compared. Finally, a brief summary is given in section 4. 2. Computational Methods 2.1. Modified BH Method. A slightly modified BH optimization method was used to obtain low-lying local minima of water clusters (n ) 20-30). In the BH algorithm,85 the energies of all randomly selected conformers are minimized before the Metropolis acceptance rule is applied, while the energy of each random conformer is directly used in the standard Monte Carlo (MC) algorithm. In fact, most of the computational cost of the BH approach is spent in the local minimization processes, where hundreds of iterations are needed at each step. However, it is not necessary to perform local minimization at each step since optimizations from neighboring configurations along the trajectory often lead to the same local minimum structure. Hence, an additional standard canonical MC (SCMC) simulation (with about 500 steps) was run to evolve the system between two neighboring local minima. In each SCMC simulation, the configuration change includes two types of trial moves for a single water molecule: (1) translation of the center of mass; (2)

Yang et al. rotation with respect to the center of mass. Each water molecule was randomly chosen with equal probability for each type of trial move. The maximum displacements of translation and rotation motion are adjustable parameters, whose values were chosen so that about 50% of the trial moves were accepted. The Euler angles were employed to describe the orientation of a water molecule in this work. In order to further accelerate the evolutive process and assist the system to cross high barriers among local minima, a higher temperature was used for the SCMC simulation. Test calculations show that 250 K is suitable for the BH process and 300 K is an appropriate temperature for the SCMC simulation. It should be mentioned that the present modification is somewhat similar to the minima-hopping method proposed by Goedecker,86 in which an additional molecular dynamics simulation was run between quenches. This modification is expected to allow the BH method to be more efficient for sampling the configurational space of water clusters. Energy minimizations were performed by using a low storage BFGS (Broyden-Fletcher-GoldfarbShanno) nonlinear optimization method.87 The convergence criterion for geometry optimizations is the root-mean-square (rms) gradient falls below 0.01 kcal/mol/Å. 2.2. AMOEBA Potential. The polarizable AMOEBA potential84 was employed in this work to describe interactions between water molecules. In comparison with nonpolarizable water potentials such as TIP3P,44 TIP4P,44 TIP5P,45 and SPC/ E,88 the major feature of the AMOEBA potential is that the fixed partial charges in traditional water models are replaced with polarizable atomic multipoles through the quadrupole moments. With the use of permanent dipoles and quadrupoles, the AMOEBA potential can reproduce accurate molecular electrostatic potentials and directional effects of hydrogen bonding and other interactions. The accuracy of the AMOEBA potential has been calibrated against high-level MP2/CBS results for some small water clusters.84 The AMOEBA binding energies are in excellent agreement with corresponding MP2/CBS results. For the hexamer clusters, the relative stabilities of several low-energy isomers computed with the AMOEBA model are also similar to those derived from MP2/CBS calculations. In addition, the geometric parameters of water clusters up to n ) 6 predicted with the AMOEBA model are generally consistent with those from MP2 calculations. These results indicate that the AMOEBA potential is suitable for studying the energies and structures of water clusters, producing more accurate results than widely used nonpolarizable water potentials. However, it should be noted that the computational cost required by this water potential is significantly higher than those needed by nonpolarizable water potentials. For example, 1000 evaluations of energy and its gradient for a (H2O)26 system (with the Tinker 4.2 package on Intel Pentium 4 3.0 GHz) need 24.53 s for the AMOEBA potential, but only 0.72 and 0.77 s for the rigid TIP3P and TIP4P potentials, respectively. Thus, a calculation with the AMOEBA potential for (H2O)26 clusters is 30 times slower compared to those with the TIP3P or TIP4P potential. 2.3. GEBF Approach. In the GEBF method,71 a target system is broken into many small fragments, and then a series of subsystems are constructed from these fragments. Each subsystem is embedded in the field of background point charges generated by all atoms outside this subsystem so that the longrange electrostatic interaction and polarization effects between remote fragments can be reasonably taken into account. Then, traditional quantum chemical calculations are carried out for these “embedded” subsystems to obtain the total energies (or

Low-Lying Structures of Large Water Clusters

J. Phys. Chem. A, Vol. 114, No. 34, 2010 9255

energy derivatives) of these subsystems, from which the total energy (or energy derivatives) of the target system can be directly evaluated. Within the GEBF approach, the total energy of a large system can be approximately expressed by71 M

ETot ≈

∑

(∑ )∑ ∑ M

CkE˜k -

k

Ck - 1

k

A

QAQB RAB A 25), we found that the lowest-energy structure occurs in one BHMC run, but similar structures with the same oxygen skeleton as the putative lowest-energy structure are produced from the other four BHMC runs. Since the number of local minimum structures for (H2O)n increases rapidly with the system size, we also added additional three independent BHMC runs for water clusters with n ) 28-30. All local minima generated from all BHMC runs were saved, and these structures are ranked in order of increasing energy. (b) For each cluster size, single-point GEBF-B3LYP/631G(d,p) calculations were first performed for the lowest 1000 AMOEBA minimum structures. Then, the lowest 100 structures were selected on the basis of their GEBF-B3LYP/6-31G(d,p) energies. For these structures, single-point GEBF-B3LYP calculations were performed again but with a large 6-311++G(d,p) basis set. (c) On the basis of the GEBF-B3LYP/6-311++G(d,p) energies, the first 10 lowest-energy structures from different structural families were used as starting points for full geometry optimizations at the GEBF-B3LYP/6-311++G(d,p) level. Here a structural family represents a subset of structures with very similar oxygen atom arrangements (but different hydrogen atom arrangements). Clearly, there is no unique way of identifying such a structural family. In this work, the two structures are assumed to belong to the same family when the rms distance of their oxygen skeletons is less than 0.5 Å. For 10 GEBFB3LYP optimized structures, we further performed vibrational frequency calculations at the GEBF-B3LYP/6-311++G(d,p) level to confirm whether these structures are true local minima, and to obtain their zero-point vibrational energies (ZPVE). In addition, for these GEBF-B3LYP optimized structures, we also ran single-point GEBF-MP2 calculations with the larger 6311++G(3df,2p) basis set to provide more accurate energies. 3. Results and Discussion Based on the AMOEBA potential, the lowest-energy structures derived from the present BH optimization procedure are shown in Figure 1 for water clusters in the range n ) 20-30. The present procedure does not guarantee that these lowestenergy structures are true global minimum structures with the AMOEBA potential, but they should be good candidates for the corresponding global minima. As shown in Figure 1, the lowest-energy structure for n ) 20 is the edge-sharing pentagonal prism, which is consistent with the global minimum structure located with TIP4P,18 SPC/E,20 and TTM2-F22 empiri-

9256

J. Phys. Chem. A, Vol. 114, No. 34, 2010

Yang et al.

Figure 1. Tentative lowest-energy structures for water clusters (n ) 20-30) obtained with the AMOEBA model. The values in parentheses are the corresponding potential energies in kcal/mol. Blue spheres represent interior oxygen atoms.

cal potentials. For n ) 21, the lowest-energy structure is a dodecahedral cage with one interior molecule, similar to the global minimum structure found with the TIP5P19 and TTM2F22 potentials. A similar clathrate structure has been experimentally observed for the H+(H2O)21 cluster, which appears as a “magic number” in many spectra.92-94 However, the lowestenergy structures predicted from the AMOEBA model are different from those obtained with the TTM2-F model22 in most cases for n g 22. For instance, the TTM2-F global minimum structure for n ) 22 has a topology of edge-sharing cubes and pentaprisms. Nevertheless, with the AMOEBA potential this structure is about 0.23 kcal/mol higher than the corresponding lowest-energy structure depicted in Figure 1. In addition, the lowest-energy structures with two internal molecules in the cage are first seen at n ) 26 for AMOEBA and TIP4P25 potentials. However, the global minimum structure of this type first occurs at n ) 28 for the TTM2-F potential.22 In comparison with conventional ab initio calculations, the GEBF approach has been demonstrated to give very accurate ground-state energies and geometric parameters of pure water clusters at various ab initio levels, as described previously71-74 and in the preceding section. First, it would be interesting to compare the performance of the GEBF-B3LYP (or GEBF-MP2) method and the AMOEBA potential in predicting the relative energies of various low-lying isomers. Table 1 provides the results for the 10 lowest-lying isomers of (H2O)26. At the AMOEBA optimized structures, the GEBF-B3LYP/6311++G(d,p) relative energies show the order 26-01 < 26-02 < 26-03, whereas a reverse ordering of 26-01 > 26-02 > 26-03 can be found from AMOEBA energies. The maximum and mean absolute errors between the GEBF-B3LYP and AMOEBA

relative energies (with respect to the energy of the 26-01 isomer) are 7.33 and 3.01 kcal/mol, respectively. On the other hand, at the GEBF-B3LYP optimized structures (with the 6-311++G(d,p) basis set), single-point GEBF-MP2 relative energies (with the larger 6-311++G(3df,2p) basis set) can be used to validate the performance of the GEBF-B3LYP method and the AMOEBA potential. With the GEBF-MP2/6-311++G(3df,2p) relative energies as the reference data, the maximum and mean absolute errors of the GEBF-B3LYP results are 4.80 and 2.08 kcal/mol, respectively, while those of the AMOEBA results are 7.13 and 2.70 kcal/mol, respectively. Hence, the results show that the accuracy of the GEBF-B3LYP/6-311++G(d,p) approach is significantly higher than that of the AMOEBA potential. This qualitative trend should hold true for other water clusters. Figure 2 shows the 10 lowest-energy structures for (H2O)26 optimized with the GEBF-B3LYP/6-311++G(d,p) approach. For other cluster sizes, the corresponding 10 lowest-lying structures (from different families) and their GEBF energies are given in the Supporting Information. Frequency calculations have verified all these structures to be local minima at the same level of theory. Since the ZPVE correction plays an important role in determining the relative stabilities of water clusters,32,95-97 the ZPVE-corrected GEBF-B3LYP energies have been employed to rank these structures. The lowest-energy structure is the isomer 26-01, which is a highly symmetrical two-centered cage structure with 3 four-membered rings, 6 five-membered rings, and 5 six-membered rings on the cluster surface. The isomer 26-02 differs from 26-01 in that its bottom surface consists of two fused four-membered rings rather than one sixmembered ring in 26-01. As seen from Table 2, the structure 26-03 is only 0.29 kcal/mol lower than 26-02, but it is a distorted

Low-Lying Structures of Large Water Clusters

J. Phys. Chem. A, Vol. 114, No. 34, 2010 9257

TABLE 1: AMOEBA, GEBF-B3LYP, and GEBF-MP2 Relative Energies (kcal/mol) (with Respect to the Energies of the 26-01 Isomer) for the 10 Lowest-Lying (H2O)26 Isomersa AMOEBA optimized geometries isomer

AMOEBA

GEBF-B3LYP/ 6-311++G(d,p)

26-01b 26-02 26-03 26-04 26-05 26-06 26-07 26-08 26-09 26-10

0 -0.918 -2.665 1.013 1.158 0.113 0.389 -2.282 2.260 3.216

0 0.559 0.906 2.001 2.189 3.161 4.764 5.050 5.174 5.309

GEBF-B3LYP optimized geometries AMOEBA

GEBF-B3LYP/ 6-311++G(d,p)

GEBF-MP2/ 6-311++G(3df,2p)

0 -1.375 -3.812 0.846 0.794 -1.183 -0.065 -4.948 -0.567 2.594

0 0.856 1.601 1.987 2.583 3.864 4.745 6.979 7.054 5.702

0 -0.918 0.526 1.617 1.438 2.588 2.507 2.182 3.245 3.423

a The structures of these isomers are optimized with the AMOEBA potential and the GEBF-B3LYP/6-311++G(d,p) approach. b For the 26-01 isomer, the AMOEBA and GEBF-B3LYP energies at the AMOEBA optimized geometries are -276.460 kcal/mol and -1988.380 155 au, respectively. The AMOEBA, GEBF-B3LYP, and GEBF-MP2 energies at the GEBF-B3LYP optimized geometries are -269.671 kcal/mol, -1988.388 896 au, and -1984.747 581 au, respectively.

Figure 2. Ten lowest-lying structures of (H2O)26 (from different structural families) optimized with the GEBF approach at the B3LYP/6-311++G(d,p) level. Blue spheres represent the interior oxygen atoms.

two-centered dodecahedral cage with 10 five-membered and 2 six-membered surface rings. In any case, all 10 lowest-energy conformers for n ) 26 display a clear two-centered cage structure. However, for n < 26, no two-centered structures are observed for the 10 lowest-energy conformers of each cluster size in this study, as shown in Figures S1-S6 of the Supporting Information. Thus, the transition from clusters with one interior molecule to clusters with two interior molecules should occur at n ) 26 from the GEBF-B3LYP calculations. Based on the ZPVE-corrected GEBF-B3LYP results, the lowest-energy structures for water clusters of different sizes (n ) 20-30) are presented in Figure 3. Meanwhile, their ZPVEcorrected GEBF-B3LYP/6-311++G(d,p) energies are listed in Table 3. To analyze the relative stability of these clusters, the binding energy per water molecule, ∆E, has also been calculated for each cluster size in Table 3. Here, ∆E is defined as the

TABLE 2: ZPVE-Corrected GEBF-B3LYP and GEBF-MP2 Energies (au) for the 10 Lowest-Lying (H2O)26 Isomersa isomer

ZPVE

ZPVE-GEBF-B3LYP/ 6-311++G(d,p)

ZPVE-GEBF-MP2/ 6-311++G(3df,2p)

26-01 26-02 26-03 26-04 26-05 26-06 26-07 26-08 26-09 26-10

0.671 561 0.672 365 0.670 704 0.670 980 0.671 070 0.670 300 0.670 756 0.672 658 0.673 016 0.671 574

-1987.717 335 (1) -1987.715 167 (3) -1987.715 641 (2) -1987.714 749 (4) -1987.713 709 (5) -1987.712 438 (6) -1987.710 579 (7) -1987.705 116 (9) -1987.704 638 (10) -1987.708 235 (8)

-1984.076 020 (3) -1984.076 679 (1) -1984.076 039 (2) -1984.074 024 (5) -1984.074 220 (4) -1984.073 156 (6) -1984.072 830 (7) -1984.071 445 (8) -1984.069 394 (10) -1984.070 552 (9)

a These structures are optimized at the GEBF-B3LYP/ 6-311++G(d,p) level, and their ZPVE values are obtained at the same level of theory. For clarity, the numbers in parentheses (n) correspond to the ordering numbers (in order of increasing energy).

9258

J. Phys. Chem. A, Vol. 114, No. 34, 2010

Yang et al.

Figure 3. Tentative lowest-energy structures for water clusters with n ) 20-30 based on ZPVE-corrected GEBF-B3LYP/6-311++G(d,p) energies. Blue spheres represent the interior oxygen atoms.

TABLE 3: Total ZPVE-Corrected GEBF Energies (au) and Corresponding Binding Energies per Molecule (kcal/mol) of the Lowest-Energy Isomer for Water Clusters (H2O)n, n ) 20-30a GEBF-B3LYP/ 6-311++G(d,p) cluster size

total energy

binding energy per molecule

20 21 22 23 24 25 26 27 28 29 30

-1529.008 068 -1605.459 076 -1681.905 803 -1758.360 343 -1834.810 658 -1911.268 742 -1987.717 335 -2064.169 276 -2140.616 802 -2217.072 039 -2293.522 847

-8.260 -8.278 -8.172 -8.289 -8.285 -8.477 -8.425 -8.455 -8.383 -8.484 -8.485

GEBF-MP2/6-311++G(3df,2p)

total energy

binding energy per molecule

-1526.205 334 -1602.518 154 -1678.827 232 -1755.139 126 -1831.453 044 -1907.767 375 -1984.076 679 -2060.388 812 -2136.701 694 -2213.014 594 -2289.326 646

-8.241 -8.317 -8.280 -8.322 -8.414 -8.509 -8.476 -8.511 -8.559 -8.605 -8.631

a The ZPVE-corrected energies for the water monomer are -76.437 242 au at the B3LYP/6-311++G(d,p) level and -76.297 134 au at the MP2/6-311++G(3df,2p) level, and the monomer geometry is optimized at the B3LYP/6-311++G(d,p) level.

energy difference between the ZPVE-corrected energy per molecule in the cluster and the ZPVE-corrected energy of a free water monomer. It can be seen from Table 3 that the stability of water clusters gradually increases with the cluster size, but no magic numbers can be obviously assigned for a given cluster size. For n ) 20, the most stable structure is the fused cubes and pentagonal prisms (20-01), which is about 1.49 kcal/mol lower in energy than the edge-sharing pentagonal prism (20-02) (see Figure S1 and Table S1 in the Supporting Information for details). In contrast to the one-centered dodecahedral cage for n ) 21, an all-surface structure is observed for n ) 22, which displays a topology of edge-sharing cubes and

pentaprisms. In general, as the cluster size further increases, various cage structures with two interior molecules among the 10 low-lying conformers can be found for n g 26. There are 2, 4, and 9 two-centered cage structures for n ) 27, 28, and 30, respectively (see Figures S7, S8, and S10 in the Supporting Information). However, the lowest-energy structure obtained here for n ) 27 is a one-centered cage fused by a pentagonal prism. Another isomer with a two-centered cage structure lies above the lowest-energy structure by 3.95 kcal/mol. For n ) 28, although the lowest-energy structure also contains only one interior molecule, the second and third lowest energy isomers, which lie above the lowest-energy structure by only 0.15 and 0.30 kcal/mol, respectively, have a two-centered cage structure. For n ) 29 and n ) 30, a two-centered cagelike structure is found to be the tentative global minimum. In particular, all 10 low-lying isomers of n ) 29 contain two interior molecules (Figure S9 in the Supporting Information). Interestingly, two interior molecules are hydrogen-bonded to each other in most cases for n ) 26 (Figure 2), but they are often separated in two-centered cage structures for n > 26. In addition, the numbers and features of hydrogen bonds in the lowest-energy structures of water clusters with n ) 20-30 described above are listed in Table 4. In general, the number of H-bonds per water tends to increase with the cluster size, but its size dependence shows an oscillatory behavior. Further, one can see that the ratio of three-coordinated water molecules to all water molecules tends to decrease from n ) 20 to n ) 30, whereas that of four-coordinated water molecules gradually increases, a trend expected for the bulk water. Nevertheless, the decreasing (or increasing) trend of three-coordinated (or fourcoordinated) water molecules is not monotonous in the cluster range n ) 20-30. Furthermore, two-coordinated water molecules seldom exist in larger water clusters (except for n ) 28 and 30).

Low-Lying Structures of Large Water Clusters

J. Phys. Chem. A, Vol. 114, No. 34, 2010 9259

TABLE 4: Number of H-Bonds per Molecule and Percentage of Two-, Three-, and Four-Coordinated Water Molecules in the Lowest-Energy Structures of Water Clusters with n ) 20-30 at the GEBF-B3LYP/ 6-311++G(d,p) Level cluster size

no. of H-bonds

two-coordinated H2O molecules, %

three-coordinated H2O molecules, %

four-coordinated H2O molecules, %

20 21 22 23 24 25 26 27 28 29 30

1.600 1.619 1.727 1.652 1.667 1.680 1.654 1.667 1.679 1.690 1.700

5.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 7.1 0.0 3.3

70.0 76.2 54.5 69.6 66.7 64.0 69.2 66.7 50.0 62.1 53.3

25.0 23.8 45.5 30.4 33.3 36.0 30.8 33.3 42.9 37.9 43.3

Since the B3LYP functional does not include the dispersion interaction, we also performed single-point GEBF-MP2/6311++G(3df,2p) calculations at the GEBF-B3LYP optimized geometries described above to obtain more accurate relative energies for the 10 low-lying structures of each cluster size. The lowest-energy structures of water clusters n ) 20-30 derived from the ZPVE-corrected GEBF-MP2 results (shown in Table 3) are shown in Figure 4. The ZPVE values used here are assumed to be the same as those obtained at the GEBFB3LYP/6-311++G(d,p) level. In comparison with the corresponding GEBF-B3LYP results (see Figure 3), six different lowest-energy structures are found for n ) 20, 23, 24, 26, 28, and 30 from the GEBF-MP2 results. For example, the most stable structure for n ) 20 is the edge-sharing pentagonal prism, which is consistent with the previous MP2 calculations with up to the aug-cc-pVQZ basis set.33 In the case of (H2O)28, the most stable structure from GEBF-MP2 calculations (Figure 4)

is a one-centered cage fused by two pentaprisms, while the lowest-energy structure from GEBF-B3LYP calculations is a one-centered cage fused by one pentaprism (with 2 twocoordinated water molecules) (Figure 3). Based on the GEBFMP2 results, the former is lower in energy than the latter by 3.05 kcal/mol. In contrast, the former lies 1.11 kcal/mol above the latter, according to the GEBF-B3LYP results (Table S8 in the Supporting Information). These differences suggest that the dispersion interaction among water molecules has a significant influence on the relative energies of low-lying isomers. 4. Conclusions In this work, the combination of the AMOEBA potential and ab initio based GEBF method has been applied to locate lowestlying structures of large water clusters (n ) 20-30). A hybrid approach is adopted for this purpose. First, a large database of low-lying isomers is obtained by using a modified BH optimization method with the polarizable AMOEBA empirical potential. Then, for each cluster size some best candidates of lowest-lying structures, screened first by single-point GEBF-B3LYP calculations, are optimized at the GEBF-B3LYP/6-311++G(d,p) level. For these GEBF-B3LYP optimized structures, vibrational frequencies have been calculated to give the ZPVE corrections. Furthermore, single-point GEBF-MP2/6-311++G(3df,2p) calculations have also been conducted at the GEBF-B3LYP optimized geometries. Several conclusions can be drawn from this study. First, although the AMOEBA empirical potential can provide good candidates for low-lying structures of water clusters, the relative stabilities of isomers predicted by this empirical potential differ noticeably from those predicted by GEBF-B3LYP/6311++G(d,p) and GEBF-MP2/6-311++G(3df,2p) calculations. Second, GEBF-B3LYP calculations indicate that for water clusters in the range n ) 20-30 a transition from one-centered

Figure 4. Tentative lowest-energy structures for water clusters with n ) 20-30 based on ZPVE-corrected GEBF-MP2/6-311++G(3df,2p) energies. Blue spheres represent the interior oxygen atoms.

9260

J. Phys. Chem. A, Vol. 114, No. 34, 2010

to two-centered cage structure first occurs at n ) 26, where all 10 low-lying structural families display various cage structures containing two interior molecules. However, the lowest-energy structures located for n ) 27 and 28 do not contain two interior molecules (some two-centered cage structures have higher energies). For n ) 29 and 30, the two-centered cage structures again dominate among the low-lying isomers including the lowest-energy structure. With increasing cluster size, the number of H-bonds per water molecule in the lowest-energy structure shows a gradually increasing trend, and the proportion of fourcoordinated water molecules gradually increases, as expected for large water clusters. Third, if the lowest-energy structures derived from GEBF-MP2 energies are compared with those from GEBF-B3LYP energies, different lowest-energy structures can be found for six cluster sizes in the range n ) 20-30, suggesting the significance of the dispersion interaction in determining the relative energies of low-lying water clusters. Thus, some recently developed DFT functionals98-101 which include the dispersion interaction may be used in the future in searching for lowestlying structures of large water clusters. On the other hand, because GEBF calculations on different subsystems are almost independent and a highly efficient parallel computation can be implemented for GEBF calculations, it will be feasible in the near future to employ ab initio based GEBF calculations with standard global optimization methods to directly locate the lowest-lying structures of large water clusters. Acknowledgment. This work was supported by the National Natural Science Foundation of China (Grants 20625309 and 20833003). Supporting Information Available: For water clusters (H2O)n for n ) 20-30 except for n ) 26, the 10 lowest-energy structures for each cluster size are given in Figures S1-S10. Additionally, the corresponding ZPVE-corrected GEBF-B3LYP and GEBF-MP2 energies are listed in Tables S1-S10. The dependence of GEBF-B3LYP results on different parameters (η) is shown in Table S11. Comparisons between GEBF-B3LYP (or MP2) and conventional B3LYP (or MP2) calculations for several water clusters are presented in Table S12. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Ludwig, R. Angew. Chem., Int. Ed. 2001, 40, 1808. (2) Hartke, B. Angew. Chem., Int. Ed. 2002, 41, 1468. (3) Nauta, K.; Miller, R. E. Science 2000, 287, 293. (4) Buck, U.; Huisken, F. Chem. ReV. 2000, 100, 3863. (5) Keutsch, F. N.; Saykally, R. J. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 10533. (6) Devlin, J. P.; Joyce, C.; Buch, V. J. Phys. Chem. A 2000, 104, 1974. (7) Devlin, J. P.; Sadlej, J.; Buch, V. J. Phys. Chem. A 2001, 105, 974. (8) Fajardo, M. E.; Tam, S. J. Chem. Phys. 2001, 115, 6807. (9) Andersson, P.; Steinbach, C.; Buck, U. Eur. Phys. J. D 2003, 24, 53. (10) Steinbach, C.; Andersson, P.; Kazimirski, J. K.; Buck, U.; Buch, V.; Ben, T. A. J. Phys. Chem. A 2004, 108, 6165. (11) Balau, L.; Wilson, K. R.; Leone, S. R.; Ahmed, M. J. Phys. Chem. A 2007, 111, 10075. (12) Mizuse, K.; Hamashima, T.; Fujii, A. J. Phys. Chem. A 2009, 113, 12134. (13) Prell, J. S.; Williams, E. R. J. Am. Chem. Soc. 2009, 131, 4110. (14) Kuyanov-Prozument, K.; Choi, M. Y.; Vilesov, A. F. J. Chem. Phys. 2010, 132, 014304. (15) Niesse, J. A.; Mayne, H. R. J. Comput. Chem. 1997, 18, 1233. (16) Guimaraes, F. F.; Belchior, J. C.; Johnston, R. L.; Robert, C. J. Chem. Phys. 2002, 116, 8327. (17) Kazimirski, J. K.; Buch, V. J. Phys. Chem. A 2003, 107, 9762. (18) Wales, D. J.; Hodges, M. P. Chem. Phys. Lett. 1998, 286, 65.

Yang et al. (19) James, T.; Wales, D. J.; Hernandez-Rojas, J. Chem. Phys. Lett. 2005, 415, 302. (20) Kabrede, H.; Hentschke, R. J. Phys. Chem. B 2003, 107, 3914. (21) Hartke, B. Phys. Chem. Chem. Phys. 2003, 5, 275. (22) Bandow, B.; Hartke, B. J. Phys. Chem. A 2006, 110, 5809. (23) Kabrede, H. Chem. Phys. Lett. 2006, 430, 336. (24) Takeuchi, H. J. Chem. Inf. Model. 2008, 48, 2226. (25) Kazachenko, S.; Thakkar, A. J. Chem. Phys. Lett. 2009, 476, 120. (26) Buck, U.; Ettischer, I.; Melzer, M.; Buch, V.; Sadlej, J. Phys. ReV. Lett. 1998, 80, 2578. (27) Tissandier, M. D.; Singer, S. J.; Coe, J. V. J. Phys. Chem. A 2000, 104, 752. (28) Lee, H. M.; Suh, S. B.; Lee, J. Y.; Tarakeshwar, P.; Kim, K. S. J. Chem. Phys. 2000, 112, 9759. (29) Lee, H. M.; Suh, S. B.; Kim, K. S. J. Chem. Phys. 2001, 114, 10749. (30) Upadhyay, D. M.; Shukla, M. K.; Mishra, P. C. Int. J. Quantum Chem. 2001, 81, 90. (31) Day, P. N.; Pachter, R.; Gordon, M. S.; Merrill, G. N. J. Chem. Phys. 2000, 112, 2063. (32) Maheshwary, S.; Patel, N.; Sathyamurthy, N.; Kulkarni, A. D.; Gadre, S. R. J. Phys. Chem. A 2001, 105, 10525. (33) Fanourgakis, G. S.; Apra, E.; Xantheas, S. S. J. Chem. Phys. 2004, 121, 2655. (34) Fanourgakis, G. S.; Apra, E.; de Jong, W. A.; Xantheas, S. S. J. Chem. Phys. 2005, 122, 134304. (35) Lenz, A.; Ojamae, L. Phys. Chem. Chem. Phys. 2005, 7, 1905. (36) Lenz, A.; Ojamae, L. Chem. Phys. Lett. 2006, 418, 361. (37) Lenz, A.; Ojamae, L. J. Phys. Chem. A 2006, 110, 13388. (38) Bulusu, S.; Yoo, S.; Apra, E.; Xantheas, S.; Zeng, X. C. J. Phys. Chem. A 2006, 110, 11781. (39) Maeda, S.; Ohno, K. J. Phys. Chem. A 2007, 111, 4527. (40) Dahlke, E. E.; Olson, R. M.; Leverentz, H. R.; Truhlar, D. G. J. Phys. Chem. A 2008, 112, 3976. (41) Nguyen, Q. C.; Ong, Y. S.; Soh, H.; Kuo, J.-L. J. Phys. Chem. A 2008, 112, 6257. (42) Hammond, J. R.; Govind, N.; Kowalski, K.; Autschbach, J.; Xantheas, S. S. J. Chem. Phys. 2009, 131, 214103. (43) Bates, D. M.; Tschumper, G. S. J. Phys. Chem. A 2009, 113, 3555. (44) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (45) Mahoney, M. W.; Jorgensen, W. L. J. Chem. Phys. 2000, 112, 8910. (46) Azhary, A. E.; Rauhut, G.; Pulay, P.; Werner, H.-J. J. Chem. Phys. 1998, 108, 5185. (47) Burant, J. C.; Strain, M. C.; Scuseria, G. E.; Frisch, M. J. Chem. Phys. Lett. 1996, 248, 43. (48) Larsen, H.; Helgaker, T.; Olsen, J.; Jørgensen, P. J. Chem. Phys. 2001, 115, 10344. (49) Niklasson, A. M. N.; Challacombe, M. Phys. ReV. Lett. 2004, 92, 193001. (50) Ochsenfeld, C.; Head-Gordon, M. Chem. Phys. Lett. 1997, 270, 399. (51) Rauhut, G.; Werner, H.-J. Phys. Chem. Chem. Phys. 2001, 3, 4853. (52) Saebø, S.; Baker, J.; Wolinski, K.; Pulay, P. J. Chem. Phys. 2004, 120, 11423. (53) Schu¨tz, M.; Werner, H.-J.; Lindh, R.; Manby, F. R. J. Chem. Phys. 2004, 121, 737. (54) Weber, V.; Niklasson, A. M. N.; Challacombe, M. Phys. ReV. Lett. 2004, 92, 193002. (55) Yang, W. Phys. ReV. Lett. 1991, 66, 1438. (56) Yang, W.; Lee, T.-S. J. Chem. Phys. 1995, 103, 5674. (57) Exner, T. E.; Mezey, P. G. J. Phys. Chem. A 2004, 108, 4301. (58) He, X.; Zhang, J. Z. H. J. Chem. Phys. 2005, 122, 031103. (59) Chen, X.; Zhang, Y.; Zhang, J. Z. H. J. Chem. Phys. 2005, 122, 184105. (60) Chen, X.; Zhang, J. Z. H. J. Chem. Phys. 2006, 125, 044903. (61) Li, W.; Li, S. J. Chem. Phys. 2005, 122, 194109. (62) Gu, F. L.; Aoki, Y.; Korchowiec, J.; Imamura, A.; Kirtman, B. J. Chem. Phys. 2004, 121, 10385. (63) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2004, 120, 6832. (64) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2005, 123, 134103. (65) Fedorov, D. G.; Ishida, T.; Uebayasi, M.; Kitaura, K. J. Phys. Chem. A 2007, 111, 2722. (66) Hirata, S.; Valiev, M.; Dupuis, M.; Xantheas, S. S.; Sugiki, S.; Sekino, H. Mol. Phys. 2005, 103, 2255. (67) Morita, S.; Sakai, S. J. Comput. Chem. 2001, 22, 1107. (68) Sakai, S.; Morita, S. J. Phys. Chem. A 2005, 109, 8424. (69) Li, W.; Li, S. J. Chem. Phys. 2004, 121, 6649. (70) Li, S.; Li, W.; Fang, T. J. Am. Chem. Soc. 2005, 127, 7215. (71) Li, W.; Li, S.; Jiang, Y. J. Phys. Chem. A 2007, 111, 2193. (72) Li, S.; Li, W. Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 2008, 104, 256.

Low-Lying Structures of Large Water Clusters (73) Li, W.; Dong, H.; Li, S. Progress in Theoretical Chemistry and Physics; Springer-Verlag: Berlin and Heidelberg, Germany, 2008; Vol. 18, p 289. (74) Hua, W.; Fang, T.; Li, W.; Yu, J.; Li, S. J. Phys. Chem. A 2008, 112, 10864. (75) Dong, H.; Hua, S.; Li, S. J. Phys. Chem. A 2009, 113, 1335. (76) Deev, V.; Collins, M. A. J. Chem. Phys. 2005, 122, 154102. (77) Collins, M. A.; Deev, V. A. J. Chem. Phys. 2006, 125, 104104. (78) Mullin, J. M.; Roskop, L. B.; Pruitt, S. R.; Collins, M. A. J. Phys. Chem. A 2009, 113, 10040. (79) Bettens, R. P. A.; Lee, A. M. J. Phys. Chem. A 2006, 110, 8777. (80) Le, H. A.; Lee, A. M.; Bettens, R. P. A. J. Phys. Chem. A 2009, 113, 10527. (81) Ganesh, V.; Dongare, R. K.; Balanarayan, P.; Gadre, S. R. J. Chem. Phys. 2006, 125, 104109. (82) Gadre, S. R.; Ganesh, V. J. Theor. Comput. Chem 2006, 5, 835. (83) Kavathekar, R.; Khire, S.; Ganesh, V.; Rahalkar, A. P.; Gadre, S. R. J. Comput. Chem. 2009, 30, 1167. (84) Ren, P.; Ponder, J. W. J. Phys. Chem. B 2003, 107, 5933. (85) Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A 1997, 101, 5111. (86) Goedecker, S. J. Chem. Phys. 2004, 120, 9911. (87) Liu, D. C.; Nocedal, J. Math. Program. 1989, 45, 503. (88) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (89) Hua, S.; Hua, W.; Li, S. J. Phys. Chem. A [Online early access]. DOI: 10.1021/jp103074f. Published Online: July 19, 2010. (90) Li, S.; Li, W.; Fang, T.; Ma, J.; Hua, W.; Hua, S.; Jiang, Y. LowScaling Quantum Chemistry (LSQC), version 2.0; Nanjing University: Nanjing, 2010. (91) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.;

J. Phys. Chem. A, Vol. 114, No. 34, 2010 9261 Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision D.01; Gaussian, Inc.: Wallingford, CT, 2004. (92) Lee, S.-W.; Freivogel, P.; Schindler, T.; Beauchamp, J. L. J. Am. Chem. Soc. 1998, 120, 11758. (93) 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. Science 2004, 304, 1137. (94) Miyazaki, M.; Fujii, A.; Ebata, T.; Mikami, N. Science 2004, 304, 1134. (95) Losada, M.; Leutwyler, S. J. Chem. Phys. 2002, 117, 2003. (96) Diri, K.; Myshakin, E. M.; Jordan, K. D. J. Phys. Chem. A 2005, 109, 4005. (97) Dahlke, E. E.; Olson, R. M.; Leverentz, H. R.; Truhlar, D. G. J. Phys. Chem. A 2008, 112, 3976. (98) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Theory Comput. 2006, 2, 364. (99) Grimme, S. J. Comput. Chem. 2006, 27, 1787. (100) Jurecka, P.; Cerny, J.; Hobza, P.; Salahub, D. R. J. Comput. Chem. 2007, 28, 555. (101) Chai, J.-D.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2008, 10, 6615.

JP1038267