5208
J. Phys. Chem. 1993,97, 5208-5210
Theoretical Study of Small Water Clusters: Low-Energy Fused Cubic Structures for (HzO)", n = 8, 12, 16, and 20 C. J. Tsai and K. D. Jordan' Department of Chemistry, University of Pittsburgh, Pittsburgh, Pennsylvania I5260 Received: January 29, 1993; In Final Form: March 30, 1993
Force field calculations employing the TIP4P potential are used to optimize several structures of the (HzO)~, n = 8, 12, 16, and 20 clusters. (H20)8is found to have, in agreement with previous studies, low-energy cubic structures of D2d and S4symmetry. The lowest-energy forms of the larger clusters are found to have fused cubic structures rather than cagelike or networked structures. The stability of the fused cubic structures of (H20)12 is confirmed by means of MP2 calculations.
The most stable form of ice, ice Ih, has a hexagonal structure with each oxygen atom surrounded by a tetrahedron of other oxygen atoms.' This structure attains the maximal degree of hydrogen bonding while being relatively free of strain. On the other hand, in finite water clusters, the maximization of the number of hydrogen bonds comes at the expenseof forming rings containing three, four, or five oxygen atoms (and the associated hydrogen atoms), introducing some strain. Thus, it is not clear that the most stable structures for small water clusters will be those that maximize the number of hydrogen bonds. It is of considerable interest, therefore, that (H20)8 in its energetically most stable form has a cubiclike structure,2d which maximizes the number of hydrogen bonds. The cubic structure reported in the theoretical studies of Brink and Glaser2 and Kim and coworkers3 has D2d symmetry. Recently, we have shown that there is a second cubic structure of S4 symmetry which is nearly isoenergetic with the D2d ~tructure.~ Figure 1 shows "ball and stick" models for these two forms of (H20),. Calculations using theTIP3Ppotentia17indicatethat these twocubic forms of (H20)8 are over 2 kcal/mol more stable than other structures5 The stability of the D2d and S4 forms of (H20), leads us to raise the question as to whether the most stable forms of (H20)12, (H20) 16, and (H20)20 might consist of "fused" cubic structures, rather than cagelike or networklike structures containing one or more five- or six-membered oxygen rings. (In this paper ring sizes are specified by the number of oxygen atoms.) To address this question, we have optimized the geometries of several forms of the (H20)8, (H20),2, (H20)16, and (H20)20 clusters. The optimizations were done with the TIP4P potential7 rather than the TIP3P potential because the former more closely reproduces the relative energies found in ab initio calculations. Initial structures, built graphically, were minimized either by means of the Monte Carlo simulated annealing procedures with sampling done using the jump-walking method9.loor by means of an eigenmode following surface walking algorithm.l0 MP2 calculations11J2 were carried out for selected TIP4P optimized structures of (H20)8and (H20) '2. These calculations employed the augmented correlation-consistent valence double-zeta + polarization (aug-ccpVDZ*) basis set.'3 (This basis set is contracted to 4s3p2d on the oxygen ,:?ms and to 3slp on the hydrogen atoms.14) The suitability of the aug-ccpVDZ* and related basis sets for describing hydrogen-bonded systems has been demonstrated in ref 15. For the (H20),cluster,we have found over 100distinct potential energy minima, in contrast to earlier studies233 which located eight or fewer minima. Figure 1 depicts, in addition to the previously discussed D2d and S4structures, the next six lowestenergy structures (in the TIP4P model), all of which are cubic, as well as one of the lowest-energynoncubic structures, designated 0022-3654/93/2097-5208$04.00/0
D2d
L1
Figure 1. Low-energy structures of (H20)8 and (H20)12obtained using the TIP4P potential.
L1. Structures with symmetry are labeled according to their symmetry point group. The TIP4P energies of these nine structures for (H20)8 and for selected low-energy forms of the (H20)12, (H20)16, and (H20)2,-~clustersaresummarizedinTables I and 11. Table I also reports MP2 energies for five of the (H20)8 structures and for three of the (H20)12structures. With the TIP4P potential, the S4structure of (H20)8 is only 0.03 kcal/mol more stable than the D2d structure but over 1.8 kcal/mol more stable than any of the other cubic structures. The MP2 calculations place the D2d structure 0.2 kcal/mol below the S4 structure and 3.9 kcal/mol below the next most stable (C2) cubic structure. These energy separations are relatively unchanged upon inclusion of correction for basis set superposition error (BSSE).I6 The energy separation between the S4 and L1 0 1993 American Chemical Society
Letters
The Journal of Physical Chemistry, Vol. 97, No. 20, 1993 5209
TABLE I: TIP4P and MPZ/aug-ccp VDZ* Energies of (H20hand (H20h cluster TIP4P (kcal/mol) MP2 (hariees)’ (H2O)s -73.027 25 -609.701 17 (-609.686 -73.053 00 -609.700 89 (-609.686 -71.205 82 -609.694 31 (-609.680 -71.178 62 -609.694 66 (-609.680 -70.664 42 -609.691 99 (-609.677 -70.632 40 -70.559 10 -70.344 60 -69.354 68 -609.691 99 (-609.677
D2d s 4
c,
c 2
CI a Cib CIc
c,
L1
____
i
81)b 51) 02) 41) 67) 4
91) r-.
(H20) I 2 -1 17.852 87 -914.565 48 (-914.544 09) -117.856 13 -1 17.859 82 -117.070 15 -1 14.609 80 -914.564 07 (-914.543 04) -114.518 34 -114.205 41 -1 14.202 64 -914.562 79 (-914.542 25)
(D2d)2 (D2d)(S4) (s4)2 (CIc) (CS)
cagel cage2 s 6
D3
The MP2calculationswerecarried out at TIP4P optimizedgeometries and excluded excitations from the 16 lowest-energyMOs of (H2O)s and 24 lowest-energy orbitals of (H~0)12. Justification for this restriction was provided by MP2 calculations on (HzO)~which showed that the relative energies of the various structures are nearly the same (agreeing on the average to 0.1 kcal/mol) when only the lowest eight M O s are frozen and when the lowest 16 M O s are frozen. bThe results in parentheses include corrections for basis set superposition error.I8
t
1.o
(1
TABLE II: TIP4P Energies (kcal/mol) of (H20)16 and (Hz0)zo (H20)16
cluster (D2d)3 (S4)(&d)(S4) (&d)(D2d)(S4) (&d)(S4)(D~d) (&d)(&)(&) (S4)3
c,
CIa CIb
network1
network2
(H20120
cluster
energy -162.861 -162.880 -162.871 -162.770 -162.777 -162.785 -161.683 -159.369 -158.876 -159.059 -159.165
78 50 83 36 42 84 54 12 63 79
(&)4 (&84)2 (s4)4
ringh
ring5b C,
network cage
energy -207.874 -207.786 -207.700 -207.350 -207.241 -206.503 -204.516 -196.789
18 30 31 94 49 95 02 25
05
structures is 3.7, 5.8, and 5.6 kcal/mol in the TIP4P, MP2, and MP2(BSSE-corrected) approximations, respectively. Eight low-energy structures for (H20)12 are also shown in Figure 1, and their energies are reported in Table 11. The three lowest-energy structures are of the fused cubic variety and are designated as (D& D284, and (S4)2,indicating their parentage in terms of the DZd and S4 forms of (H20)8. With the TIP4P potential, these three structures lie within 0.07 kcal/mol of one another, with the (S4)2 structure being most stable. There are several other fused cubic structures of (H20)12, one of which is included in Figure 1. As described by the TIP4P potential, these alternative fused cubic structures are less stable than the (D&, D2dS4,and (S4)2 structures by over 0.7 kcal/mol. The remaining low-energy structures include two cagelike structures with four five-membered rings and nine four-membered rings and two high symmetry structures with two six-membered rings and six fourmembered rings. With the TIP4Ppotentia1, these four structures are between 3.2 and 3.6 kcal/mol less stable than the lowestenergy fused cubic structure. H F and MP2 calculations with the aug-ccpVDZ* basis set were carried out for the (&)2 fused cubic structure (with overall S4 symmetry), the cagel structure, and the D3 structure. The relative energies from these and the TIP4P calculations are shown in Figure 2. The MP2/augccpVDZ* calculations place the S 4 structure lowest in energy, with the cagel and 0 3 structures lying respectively 0.9 and 1.7 kcal/mol above the S4 structure. The cagel and Ds structures
c
0
t
S4
m4P
....
MP2
w2
HF
Aug-ccpVDZ’
BSSEC
AugccpVDZ’
w BSSEC
Figure 2. Relative energies of the S4, D3, and cagel forms of (H20)12 calculated in the TIP4P model and by means of the H F and MP2 procedures using the aug-ccpVDZ*basis set and the TIP4P geometries. Ab initio results are presented with and without corrections for basis set superposition error.
are further stabilized, relative to the S4 structure, by 0.2 and 0.5 kcal/mol, respectively,upon inclusion of the correction for BSSE. Both the MP2 and TIP4P procedures predict the S4 structure to be the most stable and the Ds structure the least stable of the three structures considered. However, the energy separations between the S4 structure and the D3 and cagel structures are appreciably smaller in the MP2 approximation, indicating a deficiency in the TIP4P potential. The cagel and D3 species containing the five- and six-membered rings have two fewer hydrogen bonds than does the S4 fused cubic species. Thus, it appears that the “driving” force to maximize the number of hydrogen bonds more than compensates for the ring strain in the fused cubic forms of (H2O) I 2 . The H F procedure is inappropriate for predicting the relative stabilities of the various structures, as evidenced by the fact that it gives the opposite energy ordering for the S4, Cagel, and D3 structures than does the MP2 approximation. Figure 3 displays for (H20)16 one of the six fused cubic structures that can bederived from the DzdandS4 forms of ( H 2 0 ) 8 as well as five structures containing various numbers of fivemembered rings. Two of these, the so-called network structures, also contain one six-membered ring. The TIP4P calculations indicate that the six fused cubic structures are close in energy and also appreciably more stable than the structures containing less strained rings (but also fewer hydrogen bonds), with the energy separation between the most stable fused cubic structure and the most stable structure containing five-membered rings being about 1.2 kcal/mol. Selected optimized structures for ( H 2 0 ) 2 0 are also shown in Figure 3. These include the (D2d)4 fused cubic structures, two “ring” structures, a structure with C, symmetry, and so-called network and cage structures. The cage structure consistsentirely of five-memberedrings and has ten hydrogen atoms not involved in hydrogen bonding. The network structure has eight hydrogen atoms not involved in hydrogen bonding, whereas the fused-cubic, ring and C,structures have respectivelyfour, five, and six hydrogen atoms not involved in hydrogen bonding. With the TIP4P potential the (D2J4, (D2a4)2, and (S4)4 fused cubic structures
Letters
5210 The Journal of Physical Chemistry, Vol. 97, No. 20, 1993
We have also explored the utility of the simulated annealing method8for locating the fused cubic forms of the water clusters, starting from arbitrary structures and using the TIP4P potential. For (H20)8and (H20)12the simulated annealing method proved successfulin locating the fused cubic structures.17 However, this was not the case for (H20),6 and (H20)20, presumably because the fused cubic forms are highly disfavored on entropic grounds and, thus, were not sampled in the Monte Carlo calculations at higher tempera tures.
(D2d)3
Clb
-
Network?
Network2 -
Acknowledgment. This research was carried out with support of a grant from the National Science Foundation. The calculations were carried out on Cray YMP and Cray C90 computers at the Pittsburgh Supercomputing Center and on a large memory Cray YMP at Cray Research. We thank Carlos Sosa for carrying out the calculations on the latter systtem. References and Notes
RingSb
Network
Cage
Figure 3. Low-energy structures of (H20)1,5and (H20)200btained using the TIP4P potential.
are more stable than any of the other structures considered. The so-called network structure lies above the most stable fused cubic structure by 3.4 kcal/mol, and the cage structure is even less stable. The structures labeled by ring5a, ring5b, and Ci are respectively 0.5,0.6, and 1.4 kcal/mol less stable than the (D& cluster. Our MP2/aug-ccpVDZ* calculationson (H20) 12 indicate that the TIP4P model exaggerates the stability of the fused cubic species relative to structures containing five- and six-membered rings. It seems likely that MP2/aug-ccpVDZ* calculations on (H20)16and (H20)20 at their TIP4P geometries would locate structures with five- and six-membered rings energetically below the fused cubic structures. On the other hand, the use of MP2optimized geometries, rather than TIP4P geometries, is expected to favor the more strained fused cubic structures. Thus, although fused cubic structures may not correspond to the most stable class of structures for (H20)16 and (H20)20,they are expected to be relatively low lying in energy.
(1) Kamb, B. In Structural Chemistry and Molecular Biology; Rich, A., Davidson, N., Eds.; W. H. Freeman and Company: San Francisco, 1968; p 507. (2) Brink, G.; Glasser, L. J. Phys. Chem. 1984, 88, 3412. ( 3 ) Kim, K. S.; Dupuis, M.; Lie, G. C.; Clementi, E. Chem. Phys. Lett. 1986,131, 7766. (4) Dykstra, C. E. J. Chem. Phys. 1989,91,6472. ( 5 ) Tsai, C. J.; Jordan, K. D. J . Chem. Phys. 1991, 95, 3850. (6) Knochenmuss, R.; Leutwyler, S. J. Chem. Phys. 1992, 96, 5233. (7) Jorgensen, W. L.; Tirado-Rivers, J. J. Am. Chem. Soc. 1988, 220, 1657. (8) Kirkpatrick, S.; Gelatt, C. D., Jr.; Vecchi, M. P. Science 1983,220, 671. (9) Frantz, D. D.; Freeman, D. L.; Doll, J. D. J. Chem. Phys. 1990.93, 2769. (10) Tsai, C. J.; Jordan, K. D. Manuscript in preparation. (11) Maller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618. (12) Gaussian 92: Frisch, M. F.; Head-Gordon, M.; Trucks, G. W.; Foresman, J. B.; Schlegel, H. B.; Raghavachari, K.; Robb, M.A.; Binkley, J. S.; Gonzalez, C.; Defrees, D. J.; Fox, D. J.; Whiteside, R. A.; Seeger, R.; Melius, C. F.; Baker, J.; Martin, R. L.; Kahn, L. R.; Stewart, J. J. P.; Topiol, S.; Pople, J. A. Gaussian Inc., Pittsburgh, PA, 1991. (13) Kendall, R. A.; Dunning, T. H.; Harrison, R. J. J. Chem. Phys. 1992, 96, 6796. (14) The diffuse p hydrogen function has been deleted from the augccpVDZ* basis set, since calculations on the water dimer showed that it was relatively unimportant for the binding. (15 ) del Bene, J. In?. J. Quantum Chem., Quantum Chem. Symp. 1992, 26, 527. (16) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553; 1970, 26, 527. (17) Tsai, C. J.; Jordan, K. D. Manuscript in preparation. (1 8) For example, in estimating the BSSE for the various (HzO)xclusters, MP2 calculationswerecarried on each symmetry-unique monomer in a specific cluster both in the presence of and in the absence of the basis functions from the other centers. The BSSE for a specific water monomer is given by the differencein energies of the isolated monomer and the monomer in the presence of the basis functions on the other centers, and the total BSSE is given by the sum of the BSSE corrections for the individual monomers in the cluster.
