Clusters of glycolic acid and 16 water molecules

Chemical Physics Letters 434 (2007) 176–181 www.elsevier.com/locate/cplett

Clusters of glycolic acid and 16 water molecules Amlan K. Roy 1, James R. Hart, Ajit J. Thakkar

*

Department of Chemistry, University of New Brunswick, Fredericton, NB, Canada E3B 6E2 Received 24 October 2006; in final form 27 November 2006 Available online 8 December 2006

Abstract Semiempirical, ab initio and density functional theory calculations are used to locate many low-energy minima on the potential energy surface of the CH2OHCOOH–(H2O)16 complex. The lowest-energy structures consist of a fused, pentagonal-prism (H2O)16 cluster hydrogen-bonded to the carboxylic group of the glycolic acid; the intramolecular hydrogen bond remains intact in the glycolic acid. Low-energy structures with the acid inside a cage-like water cluster were not found.  2006 Elsevier B.V. All rights reserved.

1. Introduction Glycolic acid, CH2OHCOOH, is the simplest a-hydroxy carboxylic acid, and has both dermatological [1] and atmospheric significance [2]. Glycolic acid has rich functionality that allows it to simultaneously form intra- and intermolecular hydrogen bonds, and also allows for weaker C–H    O interactions. The rotamers of glycolic acid have been studied both experimentally [3–7] and theoretically [8,9,4,10,11,6,7]. Eight distinct rotamers have been identified. The two lowest-energy rotamers are shown in Fig. 1 as G and T. The G rotamer is the most stable and the T rotamer lies about 2.7 kcal/mol above G [11]. Both these rotamers involve an intramolecular H-bond. The most stable dimer of glycolic acid is a cyclic C2h structure consisting of two G rotamers held together by two short, linear hydrogen bonds between the carboxylic groups [12]. The next two dimer structures of higher-energy are of the same type except that one or both of the monomers are in the T rotameric form. Raman and infrared spectroscopy studies [13] suggest that glycolic acid exists in a monomeric form in dilute aqueous solution.

A first step towards understanding such a solution is to consider microsolvation in small gas-phase clusters consisting of a glycolic acid molecule (g) and a few water (w) molecules. Vibrational spectroscopy [16] has been used as a probe of water–glycolic acid interactions. Ab initio calculations [14–16] have been reported for gwn clusters with n = 1, 2, . . . , 6. The lowest-energy structures found are shown in Fig. 1. For n 6 5, they consist of a (H2O)n complex, not necessarily of lowest-energy, hydrogen-bonded to the carboxylic group of the glycolic acid. The lowest-energy gw6 structure is similar except that the water hexamer is hydrogen-bonded to both the carboxylic and a-hydroxyl groups of the acid. In all the lowest-energy clusters, the intramolecular hydrogen bond remains intact in the glycolic acid. The huge number of possible structures renders a similarly exhaustive study of significantly larger gwn clusters infeasible. We thought n = 16 was the smallest n for which cage-like structures would have a chance of being competitive with other low-energy clusters. We report a computational exploration of low-energy structures of gw16 in this Letter. 2. Computations

*

Corresponding author. Fax: +1 506 453 4981. E-mail address: [email protected] (A.J. Thakkar). 1 Present address: Department of Chemistry and Biochemistry, University of California, Los Angeles, CA 90095, USA. 0009-2614/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.12.010

As in previous work [14,15], we use a multi-stage searchand-screen method to locate local minima on the gw16 potential energy surface. In each stage, trial structures

A.K. Roy et al. / Chemical Physics Letters 434 (2007) 176–181

Fig. 1. The two lowest-energy conformers (G and T) of glycolic acid [11], and the lowest-energy structures [14,15] Gn 01 of the gwn clusters for n = 1, 2, . . . , 6.

generated by the previous stage are subject to geometry optimization at a level of theory higher than that used in the previous stage. Optimization is followed by frequency analysis to ensure that only true local minima are retained, and by single-point energy calculations in a larger basis set and/or higher-level method. Then an energy threshold is used to select a subset of the optimized structures as trial structures for the next stage. A sampling of higher-energy structures that have distinct H-bonding patterns not already present in the energy-selected set are then added to the selection. The zeroth-order step is the generation of a representative sample of guessed structures. Earlier work [14,15] on smaller gwn suggests that our search for the lowest-energy structures can be restricted to clusters in which glycolic acid is in either the G or T rotameric forms shown in Fig. 1. Binding of the water to all active sites of these rotamers was considered in a systematic fashion with greater emphasis on binding to the carboxylic group and the a hydroxyl group which are most important in low-energy structures of the smaller clusters [14,15]. Even with this strategy, the number of trial structures that need examination is extremely large. Noting that 336 guessed structures of gw6 had been used in previous work [15], we estimate that an exhaustive search for gw16 would require several thousand guessed structures. A search of that scale is not feasible with our resources. Hence, we

177

restricted ourselves to 288 guessed structures that include at least five samples of each of the structural motifs that we could think of. The 288 guessed structures were optimized using the PM3 method [17]. Single-point energies of the 285 PM3 unique local minima found were then computed at the B3LYP/6-31G**, MP2/6-31G** and B3LYP/pc1a levels. B3LYP is a hybrid DFT method [18,19] and MP2 is second-order, Møller-Plesset (MP2) [20] perturbation theory. The pc1a basis set of [4s3p1d/3s2p] contracted GTF [21,22] is a subset of Jensen’s aug-pc1 basis set [23]. Unfortunately the description in Refs. [21,22] of pc1a as having [3s1p] GTF on H atoms is incorrect. Then all 41 PM3 structures which had one or more of their three single-point relative energies less than 5 kcal/mol were selected for the second stage. Almost all of the structural motifs from the initial guessed set were present in this energy-selected subset. The main missing motif was a structure in which the acid was nearly covered by water molecules on all sides. To ensure that this and other motifs were not discarded too early, 11 additional PM3 structures were retained for the next stage. Then these 52 PM3 geometries were reoptimized at the HF/4-31G level, and three single-point energies were computed using the same methods as in the first stage. The choice of the relatively primitive PM3 and HF/431G methods for geometry optimization in the early stages was justified in earlier work on the smaller gwn complexes. We showed [14,15] that these methods lead to reasonable geometries but unreliable relative energies, and hence screening of structures should be based on single-point energies calculated using more reliable methods. As described above, this is precisely the procedure followed here. All 28 of the 52 HF local minima which had one or more of their three single-point relative energies less than 5 kcal/ mol were selected for the last stage. The ‘covered’ HF structures did not meet this energy requirement, and the lowest-energy covered structure was added to the energyselected set. These 29 HF structures were reoptimized at the B3LYP/6-31G** level. This was followed by singlepoint B3LYP/pc1a and B3LYP/pc2 energy computations where pc2 is a polarized, valence triple-zeta, [4s3p2d1f/ 3s2p1d] basis set designed specifically for DFT computations [24,25]. Harmonic B3LYP/6-31G(d,p) zero-point energy (ZPE) corrections were added to the B3LYP/pc2 energies to obtain ZPE-corrected energies denoted Ez. Similarly, ideal-gas B3LYP/6-31G(d,p) Gibbs free energy changes at 298 K and 1 atm. were shifted by the difference between the B3LYP energies in the pc2 and 6-31G(d,p) basis sets to obtain improved values of DG. Single-point B3LYP/pc2 energies, denoted Es, in the presence of a Kirkwood–Onsager reaction field with dielectric constant e = 78.39 were calculated to assess the effects of aqueous solvation. Electronic structure computations were performed with the GAUSSIAN 98 [26] and GAUSSIAN 03 [27] program suites.

178

A.K. Roy et al. / Chemical Physics Letters 434 (2007) 176–181

3. Results and discussion Table 1 lists the relative energies DE of the B3LYP/631G(d,p) local minima computed with three basis sets: 6-31G(d,p), pc1a, and pc2. Moreover, DEs includes solvation effects at the Kirkwood–Onsager level and DEz includes zero-point vibrational energy (ZPE) corrections as described above. The complexes are labeled L16 n where n increases with energy Ez. The letter L 2 {G, T, I} indicates whether the acid is in the G or T conformation, or is in an intermediate (I) conformation in which the a-hydroxyl group has rotated about the Ca–Oa bond breaking the intramolecular H-bond and possibly forming intermolecular H-bonds with water. The structures of the 15 clusters with the lowest Ez are shown in Fig. 2 and the other 14 of higher Ez in Fig. 3. Previous work on smaller gwn oligomers [15] indicates that the pc2 relative energies will provide a rank ordering of the structures that is close to the basis set limit. Table 1 shows that the relative energies computed in the smaller 6-31G(d,p) and pc1a basis sets give a different and probably less reliable rank ordering. The pc2 energies indicate

that the two lowest-energy oligomers are G1601 and G1602, and neither solvation nor zero-point vibrations change this. These two complexes have the acid molecule H-bonded to a water hexadecamer with a fused, pentagonal prism (FPP) motif. The latter is consistent with expectations based on the lowest-energy w20 clusters [28,29] which have this motif although prior, less accurate work [30] had indicated that the lowest-energy w16 and w20 oligomers consist of stacked cubes. The B3LYP/pc2//6-31G(d,p) binding energy of G1601 with respect to a glycolic acid molecule and a w16 cluster at the geometries they adopt in G1601 was found to be 19.50 kcal/mol after a counterpoise [31,32] correction of 0.49 kcal/mol for basis set superposition error (BSSE) is included. The BSSE correction is only 2.5% validating the high quality of the pc2 basis set for DFT computations. At the same level of theory, the energy required to distort the acid from its equilibrium G conformation and w16 from its optimal FPP geometry to the geometries the monomers adopt in G1601 was found to be 6.35 kcal/mol, and hence the binding energy of G1601 with respect to the G and FPP conformers is 13.15 kcal/mol.

Table 1 ˚ ) of gw16 H-bonds Relative energies (kcal/mol), and numbers and mean lengths (A Cluster

DEa

DEb

DEc

DEs d

DEz e

DGf

Ngw

ÆRgwæ

Nww

ÆRwwæ

Ntot

G1601 G1602 G1603 I1604 G1605 I1606 I1607 T1608 T1609 T1610 I1611 I1612 I1613 I1614 I1615 I1616 I1617 I1618 I1619 I1620 I1621 I1622 I1623 T1624 T1625 I1626 I1627 I1628 I1629

4.33 4.29 5.90 2.12 6.72 10.03 5.43 6.26 6.38 0.00 0.71 1.35 1.47 1.24 4.21 4.35 4.31 3.96 4.75 5.33 3.35 5.06 5.53 6.44 6.55 4.55 6.24 7.08 11.05

1.05 1.45 2.18 0.00 2.38 3.27 2.65 2.89 3.56 0.86 1.37 0.24 2.06 2.03 3.44 3.58 3.76 3.65 3.59 3.63 4.02 4.13 5.41 6.19 6.31 4.55 6.70 6.68 9.03

0.00 0.42 1.28 0.64 1.66 2.58 2.14 2.01 2.64 1.41 1.88 1.83 2.48 3.10 4.87 4.09 3.99 4.48 4.00 4.10 4.64 4.68 5.44 6.07 6.18 5.80 7.30 7.46 10.86

0.00 0.66 1.41 3.18 1.06 5.10 4.16 2.18 2.97 3.71 4.35 4.24 5.03 5.16 6.66 5.87 5.81 5.90 6.44 5.84 7.07 4.89 8.14 8.51 8.76 8.08 9.29 9.33 13.42

0.00 0.42 0.85 0.97 1.11 1.37 1.95 1.98 2.66 2.82 3.30 3.39 3.99 4.44 4.94 5.02 5.13 5.18 5.34 5.44 5.68 6.44 6.66 7.09 7.28 7.49 8.57 8.81 9.86

0.19 0.70 0.91 2.71 0.88 0.00 2.37 2.14 3.05 5.01 4.48 5.91 5.68 6.41 7.22 6.52 6.23 7.20 6.70 6.53 7.17 8.54 7.01 7.16 7.45 9.80 10.70 10.59 10.44

2 2 2 5 2 5 4 2 2 3 5 5 5 5 6 5 5 4 4 4 5 4 4 3 3 4 4 4 6

1.714 1.714 1.711 1.741 1.767 1.735 1.705 1.695 1.690 1.666 1.759 1.826 1.771 1.837 1.808 1.831 1.749 1.742 1.897 1.854 1.753 1.690 1.765 1.846 1.840 1.756 1.920 1.840 1.767

26 26 26 24 26 22 24 26 26 26 25 25 25 25 24 25 25 26 26 26 25 26 25 26 26 26 26 26 22

1.774 1.772 1.773 1.787 1.770 1.758 1.786 1.774 1.772 1.789 1.802 1.782 1.800 1.789 1.782 1.793 1.815 1.806 1.789 1.797 1.803 1.812 1.789 1.785 1.786 1.797 1.785 1.792 1.758

29 29 29 29 29 27 28 29 29 30 30 30 30 30 30 30 30 30 30 30 30 30 29 30 30 30 30 30 28

a b c d e f

B3LYP/6-31G(d,p); E(T1610) = 1527.4282143Eh. B3LYP/pc1a//6-31G(d,p); E(I1604) = 527.2879444Eh. B3LYP/pc2//6-31G(d,p); E(G1601) = 1528.2046538Eh. Onsager-B3LYP/pc2//B3LYP/6-31G(d,p); Es (G1601) = 1528.2089736Eh. Ez = E(B3LYP/pc2//6-31G(d,p))+ZPE(B3LYP/6-31G(d,p)); Ez(G1601) = 1527.7201398Eh. DG = DG(B3LYP/6-31G(d,p)) + E(B3LYP/pc2//6-31G(d,p))E(B3LYP/6-31G(d,p)); relative to DG(I1606) = 1527.7938570Eh.

A.K. Roy et al. / Chemical Physics Letters 434 (2007) 176–181

C

179

C

C

C

C

C

I1615

I1614

I1613

C

I1612

C

C C

C

I1611 C

T1610

C

C C

C

C

C

T1608

T1609

I1607

c

c

C

C

I1604

C

G1601

C

G1605

C

C

C

C

G1602

I1606

C

C

G1603

Fig. 2. The most stable gw16 clusters numbered in order of increasing energies including the zero-point correction.

Fig. 1 summarizes previous work [15] which showed that only the carboxylic group of the acid is H-bonded to a wn cluster in the lowest-energy structures of gwn with n 6 5 but both the carboxylic and a-hydroxyl groups are H-bonded to a water hexamer in the lowest-energy gw6 complex. Although we had expected the lowest-energy gw16 clusters to be like gw6 in this respect, G1601, G1602 and G1603 have only the carboxylic group of the acid H-bonded to a w16 complex. However, our expectations are met when the Gibbs free energy is considered, and I1606 is the most favored structure. It involves a water decamer in the form

of a pentagonal prism attached to a water hexamer not unlike the book form. The acidic hydroxyl H-bonds to the prism whereas Oa and the a-hydroxyl group H-bond to the hexamer. This structure has only 27 H-bonds of which only 22 are between pairs of water molecules. The intramolecular H-bond in the acid moiety remains intact in nine structures, including seven of the 10 lowestenergy ones, but notably not in I1606. Table 1 lists the total number of intra- and intermolecular H-bonds (Ntot), the number of H-bonds between water molecules (Nww), and the number of H-bonds between the acid and water

180

A.K. Roy et al. / Chemical Physics Letters 434 (2007) 176–181

c

I1629

c

C

C

I1628

I1626 C

C

T1625 C

I1627

C

C

C

I1622

T1624

I1623

C C C

C

C C

C

C

C

C

I1620

C

C

I1621 I1619

I1616

C

C

C

C

I1617 C

C

I1618

Fig. 3. Other gw16 complexes numbered in order of increasing energies including the zero-point correction.

molecules (Ngw). Two-thirds (18) of the clusters have Ntot = 30 and eight others have Ntot = 29, whereas three clusters (I1606, I1607 and I1629) have fewer H-bonds.

All but five gw16 complexes have either 25 or 26 water– water H-bonds. The number of H-bonds between the acid and water molecules shows the greatest variability with Ngw

A.K. Roy et al. / Chemical Physics Letters 434 (2007) 176–181

ranging from 2 to 6. The nine lowest-energy structures all have Ntot < 30 showing that the structures with the most H-bonds are not necessarily the most stable ones. The relative stability of a structure clearly depends on the strength as well as the number of H-bonds. An indication of the average strength of the H-bonds is given in Table 1 by the mean H . . . O distances ÆRgwæ and ÆRwwæ for acid–water and water–water H-bonds, respectively. Quantitative correlations between relative energies and simple measures such as Nww, Ngw, ÆRgwæ and ÆRwwæ are not easy to find. Next, we consider some of the structural motifs in the other gw16 clusters. Two different types of FPP structures are seen. One can be thought of as consisting of layered, bicyclic octagons of water molecules; see, for example, G1601 and G1602. The other motif, seen in G1603 and G1605, has one of the pentagonal prisms rotated by roughly 90relative to the other prism. A pentagonal prism-shaped decamer and a hexamer bonded to the carboxylic and a-hydroxyl groups respectively are seen in several complexes including I1604, I1606, and I1607. Two or more fused cubes of water molecules are wide-spread in the higher-energy oligomers including I1610, I1613, I1620 and I1627. A cube fused to a pentagonal prism is seen in I1619, for example. Cage-like structures with the acid inside are not energetically competitive. After optimization at the B3LYP/6-31G(d,p) level, the most covered HF/431G structure we found turned into I1629 where the acid has transferred a proton to the water molecules. Finally, our search-and-screen method was tested by raising the first-stage energy selection threshold from 5.0 kcal/ mol to 5.3 kcal/mol. The nine additional PM3 structures that met this less stringent criterion were carried through all remaining stages without any further screening. The final DEz values of these structures ranged between 2.96 kcal/mol and 11.4 kcal/mol. Comparison with Table 1 shows that the search-and-screen strategy with the original 5 kcal/mol threshold correctly identified the ten lowest-energy structures. The additional structures all contained two fused cubes. 4. Final remarks A search-and-screen strategy employing semiempirical, ab initio, and density functional methods was used for a preliminary, non-exhaustive exploration of the potential energy surface of gw16. Low-energy structures with the acid inside a cage-like water cluster were not found. The lowestenergy structures consist of a fused, pentagonal-prism (H2O)16 cluster hydrogen-bonded to the carboxylic group of the glycolic acid which retains its intramolecular hydrogen bond. A structure in which the a-hydroxyl and carboxylic groups of the acid bind to different parts of a water cluster had the lowest Gibbs free energy. Definitive determination of the lowest-energy structure will require the use of high-level ab initio methods.

181

Will gwn clusters with even more water molecules be analogous? Will their lowest-energy structures consist of a glycolic acid molecule bound to the exterior of a wn oligomer? A study of gw28 clusters is being carried out in our laboratory to address these questions. Acknowledgements We thank Catherine Nordstrom for help with figure preparation, and Shaowen Hu for preliminary discussions on possible structures of the slightly smaller gw14 clusters. This work was supported by the Natural Sciences and Engineering Research Council of Canada. References [1] A.D. Katsambas, Dermatology 210 (2005) 39. [2] S.R. Souza, P.C. Vasconcellos, L.R.F. Carvalho, Atmos. Environ. 33 (1999) 2563. [3] C.E. Blom, A. Bauder, J. Am. Chem. Soc. 104 (1982) 2993. [4] H. Hollenstein, T.-K. Ha, H.H. Gu¨nthard, J. Mol. Struct. 146 (1986) 289. [5] K. Iijimaa, M. Katoa, B. Beagley, J. Mol. Struct. 295 (1993) 289. [6] P.D. Godfrey, F.M. Rodgers, R.D. Brown, J. Am. Chem. Soc. 119 (1997) 2232. [7] I.D. Reva, S. Jarmelo, L. Lapinski, R. Fausto, Chem. Phys. Lett. 389 (2004) 68. [8] M.D. Newton, G.A. Jeffrey, J. Am. Chem. Soc. 99 (1977) 2413. [9] T.-K. Ha, C.E. Blom, H.H. Gu¨nthard, J. Mol. Struct. (Theochem) 85 (1981) 285. [10] M. Flock, M. Ramek, Int. J. Quantum Chem. S26 (1992) 505. [11] F. Jensen, Acta Chem. Scand. 51 (1997) 439. [12] N.E.-B. Kassimi, E.F. Archibong, A.J. Thakkar, J. Mol. Struct. (Theochem) 591 (2002) 189. [13] G. Cassanas, M. Morssli, E. Fabreque, L. Bardet, J. Raman Spectrosc. 22 (1991) 11. [14] A.J. Thakkar, N.E.-B. Kassimi, S. Hu, Chem. Phys. Lett. 387 (2004) 142. [15] A.K. Roy, S. Hu, A.J. Thakkar, J. Chem. Phys. 122 (2005) 074313. [16] D.K. Havey, K.J. Feierabend, K. Takahashi, R.T. Skodje, V. Vaida, J. Phys. Chem. A 110 (2006) 6439. [17] J.J.P. Stewart, J. Comput. Chem. 10 (1989) 209. [18] A.D. Becke, J. Chem. Phys. 98 (1993) 5648. [19] P.J. Stephens, F. Devlin, C. Chabalowski, M. Frisch, J. Phys. Chem. 98 (1994) 11623. [20] Chr. Møller, M.S. Plesset, Phys. Rev. 46 (1934) 618. [21] A.K. Roy, A.J. Thakkar, Chem. Phys. 312 (2005) 119. [22] J.R. Hart, A.J. Thakkar, J. Mol. Struct. (Theochem) 714 (2005) 217. [23] F. Jensen, J. Chem. Phys. 117 (2002) 9234. [24] F. Jensen, J. Chem. Phys. 116 (2002) 7372. [25] F. Jensen, J. Chem. Phys. 118 (2003) 2459. [26] M.J. Frisch et al., GAUSSIAN 98, Revision A.11, Gaussian Inc., Pittsburgh, PA, 2001. [27] M.J. Frisch et al., GAUSSIAN 03, Revision B.05, Gaussian Inc., Pittsburgh PA, 2003. [28] G.S. Fanourgakis, E. Apra`, S.S. Xantheas, J. Chem. Phys. 121 (2004) 2655. [29] A. Lenz, L. Ojama¨e, Chem. Phys. Lett. 418 (2006) 361. [30] S. Maheshwary, N. Patel, N. Sathyamurthy, A.D. Kulkarni, S.R. Gadre, J. Phys. Chem. A 105 (2001) 10525. [31] S.F. Boys, F. Bernardi, Mol. Phys. 19 (1970) 553. [32] I. Mayer, P.R. Surja´n, Chem. Phys. Lett. 191 (1992) 497.