Two-Dimensional Icy Water Clusters Between a Pair of Graphene-Like

Aug 3, 2016 - In this case, a periodic cell of graphene substrate must be large enough to avoid ...... Binding Energies (BEs, kcal/mol) of (H2O)1–6 ...
0 downloads 0 Views 5MB Size
Download PDF
Article pubs.acs.org/JPCC

Two-Dimensional Icy Water Clusters Between a Pair of GrapheneLike Molecules or Graphene Sheets Seong Kyu Kim,*,† Wenzhou Chen,‡ Saeed Pourasad,‡ and Kwang S. Kim*,‡ †

Department of Chemistry, Sungkyunkwan University, Suwon 16419, Korea Center for Superfunctional Materials, Department of Chemistry, Ulsan Institute of Science and Technology (UNIST), Ulsan 44919, Korea

‡

S Supporting Information *

ABSTRACT: To understand the 2-dimensional (2D) structural evolution of water molecules intercalated into a graphene bilayer, the geometries of water clusters up to tridecamer formed between a pair of graphene sheets or between graphene-like molecules (coronene and dodecabenzocoronene) are investigated. Due to their large sizes, the selfconsistent-charge density-functional tight-binding (SCCDFTB) method expanded into the third order and supplemented with a Slater−Kirkwood dispersion term was used. In this way both hydrogen bonding and H−π/π−π interactions are calculated in a balanced manner with the right magnitude of binding energies very close to the reference values based on the most accurate ab initio results. It should be noted that conventional density functional theory (DFT) calculations underestimate the H−π/π−π interaction, while dispersion-corrected DFT calculations overestimate hydrogen bonding. The latter method is also employed for comparison and to confirm the reliability of the SCC-DFTB results. For (H2O)6, the fused bitetragonal hexamer is nearly isoenergetic to the most stable planar hexagonal ring structure, and it is more frequently found. In (H2O)10 and (H2O)13 clusters, a tetragon is the most frequent geometry followed by a pentagon, while the hexagon is less frequent. These results certainly provide evidence of the recent planar tetragonal ice structure found inside a graphene bilayer (Nature 2015, 519, 443) which is in contrast to the well-known hexagonal pattern of the bulk ice. This structural change from hexagonal to tetragonal network on the graphene surface is attributed mainly to the inherent nature of 2D water.

1. INTRODUCTION Two-dimensional (2D) chemistry has been emerging through noncovalent adsorption on graphene1 or intercalation inside bilayer graphene.2 It is thus timely to study water structures in 2D systems since water is the most important substance in our life. As yet, such the 2D nature of water has rarely been studied. Very recently, it has been reported that ice water forms a 2D checkerboard pattern between a pair of graphene sheets.3 The checkerboard arrangement of water molecules should withstand substantial strain out of unusual hydrogen bonding structure. The understanding of such complex phenomena cannot be made simply either by ab initio or DFT calculations or by empirical/ semiempirical potential-based Monte Carlo or molecular dynamic (MD) simulations. For the logical understanding, it is necessary to study the progressive change as the number of water molecules increases. Icy water confined in a highly restricted environment shows unique properties that are not observed in bulk ice. The confined water molecules can assemble into structures that are different from their tetrahedral hydrogen bonding network of liquid or ice. When the confined environment is hydrophobic, they usually follow the confined shape. For example, according to MD simulations4 and X-ray diffraction studies,5 icy water forms tubular structures inside carbon nanotubes with tetragonal to octagonal cross sections depending on the tube diameter. When © XXXX American Chemical Society

confined by 2D substrates, a monolayer ice is likely to form a 2D structure.6 As already mentioned, icy water can form a 2D checkerboard pattern between a pair of graphene sheets.3 It is intriguing because water tends to favor hexagon structure in ice and liquid, and the pentagon structure is relatively stable in the gas phase, while the planar tetragon structure is strained. In this regard, to understand the 2D icy water structure, a systematic study for water−water and water−substrate interactions is required with accurate calculations. However, high level calculations are limited to small systems. Calculations of large systems are feasible only at low levels of theory. Here, the information from small systems may help to set up initial conditions for large systems. With these in mind, we have tried to find the structures and interactions of water clusters between a pair of graphene-like polyaromatic rings or a pair of graphene sheets. Two kinds of approaches are used to find the information on water−graphene interactions. The first is to use graphene directly as a substrate. In this case, a periodic cell of graphene substrate must be large enough to avoid interaction of water molecules between neighboring cells. This involves a large number of atoms Received: June 16, 2016 Revised: July 31, 2016

A

DOI: 10.1021/acs.jpcc.6b06076 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Table 1. Binding Energies (kcal/mol) of Selected Clusters Calculated at Various Levels of Theorya SCC-DFTB (H2O)2 (H2O)3 (H2O)4 (H2O)5 (H2O)6, ring (H2O)6, book (H2O)6, bag (H2O)6, cage (H2O)6, prism Bz(H2O)1 Bz−Bzb Cor−Corb Cor−Cor (6 Å)c

PBE-D3

3D

2D

3

TZVP (scaledd)

TZVPP (scalede)

ref values

5.05 15.8 27.5 36.3 44.8 45.9 46.1 47.0 48.0 2.90 3.17 19.3 1.55

3.40 9.84 18.2 23.8 29.1 30.2 30.0 30.6 31.0 1.79 3.21 20.1 0.75

4.97 15.6 26.9 35.5 43.8 44.8 45.1 45.9 46.6 2.17 0.19 0.0 0.0

7.02 (4.91) 22.3 (15.6) 39.1 (27.4) 51.4 (36.0) 63.2 (44.2) 65.1 (45.6) 64.8 (45.4) 66.1 (46.3) 66.7 (46.7) 5.19 (3.63) 2.90 (2.03) 17.6 (12.3) 1.86 (1.30)

6.31 (5.05) 20.4 (16.3) 35.6 (28.5) 46.6 (37.3) 56.9 (45.5) 59.2 (47.4) 59.0 (47.2) 60.3 (48.2) 60.6 (48.5) 4.55 (3.64) 2.88 (2.30) 16.8 (13.4) 1.93 (1.54)

5.0,f 5.02g,h 15.8f,i,j 27.4,f 27.6j 35.9,f 36.3j 44.3,f 44.9,j 45.2k 45.4,f 45.6,j 46.1k 45.5k 45.9,f,j 46.5k 46.2,f 45.9,j 46.6k 3.34g,l 2.7−2.8m 20.0 ± 7.3,n 20.0 ± 1.7o

a

All energies at the SCC-DFTB and PBE-D3 levels are given without corrections for zero point energy and basis set superposition error (BSSE). Parallel-displaced geometry. cCor−Cor (6 Å) denotes that two Cor molecules are separated by 6 Å, where the stacked or the parallel-displaced conformations do not show any significant difference in binding energy. dScaled binding energies in parentheses scaled by 0.7 toward the realistic values because of large overestimation. eScaled binding energies in parentheses scaled by 0.8. fCCSD(T)/CBS result from ref 31. gCCSD(T)/CBS result from ref 32. hCCSD(T)/CBS result from ref 33. iMP2/CBS results from ref 34. jMP2/CBS result from ref 35. kMP2(FC)/HZ4P(2fg,2d)++ result from ref 36. lCCSD(T)/CBS result from ref 37. mCCSD(T)/CBS result from refs 38−40. nSCS-MP2/pVDZ result calculated in this work (the middle value between BSSE uncorrected and corrected values, where the error bar denotes half the BSSE; see refs 28, 36, and 42 for the reliability of the middle value). oM06-2X/DIDZ result from ref 41 (the middle value between BSSE uncorrected and corrected values, where the error bar denotes half the BSSE). b

density functionals would not be feasible for this systematic study involved in very large systems. What is more critical is that dispersion-corrected DFT-D3 gives overbinding energies (as much as 20−30%) for water clusters. Therefore, here we exploit a self-consistent-charge density-functional tight-binding (SCCDFTB) method21 which provides more realistic binding energies. The SCC-DFTB method is an approximation to DFT by expanding DFT energy expression to second or third order terms. The interaction between atoms is then described by an effective damped Coulombic interaction between atomic charges derived from the expansion. Since the speed of calculation is higher by 2−3 orders of magnitude than that of DFT calculations using small to medium size basis sets and the accuracies in describing molecular geometries are comparable to the results of DFT or higher level calculations, the SCC-DFTB is highly useful for large systems. As most density functionals describe dispersion interaction poorly, additional terms with empirical fitting parameters are usually employed to reproduce the dispersion interaction. The use of empirical fitting parameters provides flexibilities to the SCC-DFTB method rather than lessens its reliabilities, as the parameters can be slightly adjustable to best describe the system under study. In our case, a Slater−Kirkwood atomic polarizable model was used with the parameters described by Elstner, Hobza, and co-workers,22 but a slightly larger (by 0.2 Å) covalent radius for C atom was used to provide better dispersion energies for the benzene (Bz) dimer and Bz−water systems. In the standard SCC-DFTB method, the second order expansion of the DFT energy is used, and this is sufficient for most covalently bound systems. However, it is insufficient for hydrogen bonding systems, which require the expansion up to the third order.23 To be short, the major computation method employed in this work is referred to as SCC-DFTB-3D, where “3” represents the expansion up to the third order and “D” represents the Slater−Kirkwood type dispersion term. These features are provided by a version 1.2 code of the DFTB+

in the calculation, and therefore, the level of calculation is rather limited. Limited studies7−10 based on density functional theory (DFT) have been performed on water clusters up to hexamer on a graphene monolayer. The second kind of approach is using polycyclic aromatic hydrocarbons as a model substrate, and then, the results from the study are extrapolated to that of graphene.11−19 Reported adsorption energies of a single water molecule on graphene fall in a wide range (between 2.2 and 5.8 kcal/ mol).11−19 Despite the large uncertainty in the reported adsorption energy, the geometry of the adsorbed water is relatively consistent from method to method. On the most stable geometry, the water molecule stands upright with the two hydrogens pointed down toward the substrate plane, although the geometry with the two hydrogens pointed up is only slightly higher in energy (by ∼0.5 kcal/mol). Compared to many computational works for water on graphene or its molecular models, the study of water between two layers of graphene or its molecular model is very rare. Ruuska and Pakkane20 used a Hartree−Fock (HF) method for a single water molecule placed between a pair of coronene (Cor; C24H12) molecules. However, the water intercalation energy could not be correctly described due to the lack of electron correlation. In this article, we report the energies and geometries of (H2O)1−6 clusters between a pair of coronene molecules, (H2O)n=1−6,10 clusters between a pair of dodecabenzocoronene (DBC; C54H18) molecules, and (H2O)n=10,13 species between a pair of graphene sheets, with a hope that such information would provide geometries of much larger water clusters inside a bilayer graphene.

2. METHODS Since accurate ab initio calculations including coupled cluster with single, double, and perturbative triple excitations (CCSD(T)) are too time-consuming, the only feasible approach would be to use DFT with dispersion correction. Even the hybrid B

DOI: 10.1021/acs.jpcc.6b06076 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C quantum mechanical calculation package.24,25 DFT calculations were carried out using the Turbomole program package (version 6.4).26 All the geometrical structures of water clusters with and without rigid substrates were optimized using both SCC-DFTB3D and the PBE-3D/TZVP methods independently. 2.1. Background Study: (H2O)1−6, H2O−Benzene, Benzene Dimer, and Coronene Dimer. The geometries of water clusters confined by a pair of graphenes or graphene-like molecules are determined mostly by H−OH H-bonding interaction in water clusters and by H−π interaction27 between water and substrate, while the π−π interaction28 between the substrate pair contributes to a lesser extent. In Table 1, the binding energies of (H2O)1−6, Bz(H2O)1, Bz dimer, and Cor dimer calculated by the SCC-DFTB methods are compared with those of the highest level ab initio calculations from the literature and those of PBE29-D330 calculations which are commonly and widely used in both physics and chemistry communities as one of the most representative dispersion-corrected DFT methods. Three values listed under columns of 3D, 2D, and 3 are, respectively, the binding energies from the SCC-DFTB calculations with both the third order expansion and the dispersion term (3D), with the second order expansion and the dispersion term (2D), and with the third order expansion term without dispersion (3). It should be noted that PBE-D3 overestimates the binding energies; thus, the energy values are scaled to obtain almost the right values calculated at the high level of ab initio theory including nearly full electron correlation. In contrast, the binding energies of water clusters calculated with SCC-DFTB-3D are almost same as the most accurate reference values within uncertainties of the reference methods. In water clusters, the contribution of the dispersion is a small fraction. In contrast, the expansion up to the third order is very important; without the third order expansion, the binding energies of water clusters are underestimated as much as about 70%. In the binding energy of Bz(H2O)1 calculated by SCC-DFTB-3D, the dispersion term contributes about 25% to the total binding energy, which is also in good agreement with the reference values. In the Bz dimer and the Cor dimer, the dispersion term is dominant in the binding energy. The binding energy of the Bz dimer calculated by the SCC-DFTB-3D is larger by ∼10% than the reference value. However, this small overestimation hardly influences the geometries of water clusters in the substrates of this work since the π−π interaction is very weak as the intersubstrate distance is far from their binding distance due to the intercalated water clusters. From these results, SCC-DFTB-3D performs very well, in close agreement with the best reference values since water− water, water−π, and π−π interactions are all in reasonable agreement with the reference values.31−41 The results are much more reliable than both the standard DFT (which gives reasonable values of water−water interactions but poor values for water−π, and π−π interactions) and the standard DFT-D3 (which gives overbinding energies for water−water interactions despite reasonable values for water−π, and π−π interactions). As noted in Table 1, the PBE-D3/TZVP and PBE-D3/TZVPP binding energies need to be scaled by 0.7 and 0.8, respectively, so as to give the realistic magnitude of these water binding energies. Furthermore, their water−water and water−Bz binding energies are much less reliable than the SCC-DFTB-3D results, when compared with the accurate CCSD(T)/complete-basis-set(CBS) results. 2.2. Preparing Initial Geometries. The three substrates employed in this study are a pair of Cor molecules, a pair of DBC

molecules, and a bilayer graphene. The structure of a single Cor molecule was fully optimized to be used as a molecular substrate. Another molecular substrate DBC is partially optimized with the central six C−C bond lengths fixed to 1.427 Å. A single layer of graphene consisting of 128 carbon atoms was constructed with all C−C bond lengths fixed to 1.427 Å in a periodic supercell of 19.773 × 17.124 × (20 + dz) Å3, where dz is the interlayer separation. The internal geometries of the three kinds of substrates were then fixed in the next step calculations for the intercalated water systems. The general approach of computation is as follows. First, each water cluster (H2O)n (n = 1−6) on the surface of a single Cor was optimized. Here, the initial geometry of each water cluster on Cor was taken from its known geometry in an isolated state.36,41,43 All initial geometries of water clusters were positioned approximately at the center of the Cor plane. Planar ring structures for n = 3−6 were positioned with the ring plane parallel to the Cor plane. Each nonplanar structure for n = 6, such as book, bag, cage, and prism configurations,43 was positioned with two or three different tumbling orientations and two lateral spinning orientations with respect to the Cor plane. Second, the other Cor was placed above each optimized (H2O)1−6/Cor geometry, and then the water cluster between the Cor pair was optimized. The two Cor molecules were stacked with the Cor− Cor distance (dz) set to 20 Å initially, and then, the distance was reduced gradually. The geometry of water cluster optimized at each dz was used as an initial geometry for optimization at a next decremented dz. As the energy of the cluster approaches a minimum, the decrement of dz was reduced to as much as 0.1 Å. The M06-2X/DIDZ41 calculations show that the paralleldisplaced Cor dimer (displaced by 1.75 Å with the interlayer distance of 3.32 Å) is more stable than the exactly stacked sandwiched Cor dimer (with the interlayer distance of 3.66 Å). However, at the separation distance of 6 Å, we find that both stacked and parallel-displaced structures have almost the same binding energies, and so the orientation is not important. A similar trend is also found for the SCC-DFTB-3D method which shows the binding energies of 19.3 kcal/mol for the paralleldisplaced Cor dimer (3.3 Å) and 1.6 kcal/mol for the Cor−Cor (6 Å) conformer. When any water clusters in this study are inserted into the Cor pair, the binding energy between two Cor molecules at dz = 6.0 Å is almost the same (within 0.07 kcal/mol) regardless of their stacking orientations. This trend also holds for the DBC pair and the bilayer graphene. Since such stacking orientations do not affect the binding energies significantly, the water structures could change the stacking conformation in favor of their maximal interactions with two Cor molecules, which, however, turn out to be insignificant. Overall, the energy difference depending on the stacking types could be minimal. Nevertheless, we calculated the stacking cases and displaced stacking cases, and the intermediate conformations between the two are also studied; the maximum binding energies are reported here. Third, a pair of DBC molecules were used as the substrate to adopt a larger size (H2O)10 as well as (H2O)1−6. (H2O)10 was chosen because a bihexagonal configuration, like naphthalene, can be formed. As configurations of various water clusters are diverse, we will use the following notation. The ring type is denoted as 3 for triangle, t for tetragon, p for pentagon, and h for hexagon. The cross(es) ( × ) separates the ring types between successive rows. When similar configurations are to be distinguished, additional italic letters, a, b, and so forth, are attached at the end. We will use the notation throughout for all C

DOI: 10.1021/acs.jpcc.6b06076 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

3. RESULTS 3.1. (H2O)1−6 on Graphene-like Surface. In our optimized geometry of (H2O)1 on graphene-like surface, the two hydrogens point down (d) toward the substrate plane, consistent with most results in the literature; this geometry will be denoted as “dd”. The geometries of (H2O)2−6 on the Cor surface optimized with the SCC-DFTB-3D are shown in Figure 2. The geometries on

fused ring configurations in this study. The configurations and the notations of (H2O)3−6,10,13 clusters are schematically represented in Figure 1.

Figure 1. Schematic representation of the hydrogen bonding network configurations of 2D water clusters in this study.

An initial geometry of the bihexagonal (hh) configuration was created under the following guideline; 11 hydrogen bonds are needed to form a bihexagon, and at a stable geometry the remaining 4 and 5 OH bonds point alternatively toward the top and bottom DBC molecules, respectively. Four checkerboard configurations that consisted of four fused tetragonal rings are possible for (H2O)10. They are ttt × t-a, ttt × t-b, tt × tt, tttt, as shown in Figure 1. The following guideline is applied to prepare initial geometries for the checkerboard configurations; 13 hydrogen bonds are needed to keep the structure, and the remaining 3 and 4 OH bonds are arranged to point toward the top and bottom DBC molecules, respectively. Configurations can also be created with combinations of fused triangular, tetragonal, and pentagonal rings. Initial geometries of such kinds of structures were found with a simulated annealing method.44 These structures (tt × tp, pt × p, tpp, tht-a. tht-b, h × tp) are depicted in Figure 1. A series of MD simulations were carried out starting at 400 K and ending at 0 K in order to overcome conformational barriers followed by gradual cooling to reach a low energy regime. Annealing was done with Nose− Hoover thermostat and 0.5 fs time step. During the simulation, we used the TIP4P model for water. A pair of DBC molecules with the interlayer separation varied between 6.0 and 7.0 Å was used as a substrate to find six stable configurations of (H2O)10. Fourth, to adopt (H2O)13 clusters, as well as (H2O)10 clusters mentioned above, a bilayer graphene was used as a substrate. This size of the water clusters was chosen because a trihexagonal configuration (hh × h) can be formed. Seven checkerboard configurations, which consist of six fused tetragonal rings, were also prepared for (H2O)13. Configurations consisting of a combination of multiple kinds of rings are also possible. Initial geometries in such configurations were prepared by the simulated annealing method.

Figure 2. Top view of Cor(H2O)2−6 optimized at the SCC-DFTB-3D level. For Cor(H2O)6, only the structures with the lowest energy at each nonplanar configuration (ring, book, bag, cage, prism) are displayed. Hydrogen bonds between water molecules are shown with green dotted lines.

other graphene-like surfaces are qualitatively independent of the substrate studied in this work. The geometry of the lowest energy water dimer shows C2v point group symmetry with the principal axis along the H-bond. All three non-H-bonding hydrogens point down toward the Cor surface. One dangling OH of the H-bond donating water molecule makes 14.3° from the surface normal and the two dangling OH’s of the H-bond accepting water molecule make 62.3° from the surface normal. We call this structure “dHdd”, which is ∼0.4 kcal/mol more stable than the structure “uHdd” with the first dangling hydrogen (Hd) pointing upward (u), where “H” denotes a H-bonded hydrogen. The geometries of (H2O)3−5 take on planar ring structures, forming an equilateral triangle, a square, and a regular pentagon, respectively, consistent with the results in the literature.10,17 The geometries of water hexamer on graphene-like surface take on complicated features as the isolated water hexamer may exist in five different isoenergetic geometries: cyclic ring, book, bag, cage, and prism. Lin et al. have studied, with the SCCDFTB-D method, the geometry changes of the water hexamers as they are adsorbed onto the graphene-like surface.17 Our results recognize the following facts which are consistent with their discussion points. (i) In all kinds of water hexamers, geometry change upon adsorption is rather small. The negligible geometry change includes the orientation of dangling OH’s with respect to the skeleton of water cluster. As a result, the number of dangling OH bonds pointing toward the substrate surface is a major factor to determine relative stability of the corresponding hexamer. (ii) D

DOI: 10.1021/acs.jpcc.6b06076 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Table 2. Binding Energies and RMSD Values of (H2O)1‑6 on a Single Layer of Graphene-like Substrate (Cor, DBC, Graphene), Optimized with the SCC-DFTB-3D and PBE-D3/TZVP Methodsa binding energy (kcal/mol)a (H2O)1 (H2O)2 (H2O)3 (H2O)4 (H2O)5 (H2O)6, ringb (H2O)6, bookb (H2O)6, bagb (H2O)6, cageb (H2O)6, prismb

RMSD (Å)

Cor

DBC

graphene

Cor

DBC

graphene

3.00 [5.41 (3.79)], {2.4−5.2}c 8.29 [16.2 (11.3)] 20.3 [32.4 (22.7)] 32.8 [49.8 (34.9)] 43.5 [63.6 (44.5)] 52.8 [78.5 (55.0)] 53.8 [77.6 (54.3)] 53.1 [78.5 (55.0)] 52.8 [76.6 (53.6)] 53.3 [76.7 (53.7)]

2.75 [5.67 (3.97)], {2.6−5.3}d 9.49 [16.0 (11.2)] 22.1 [33.1 (23.2)] 35.0 [51.0 (35.7)] 45.8 [65.5 (45.9)] 55.6 [80.7 (56.5)] 56.3 [79.4 (55.6)] 55.8 [79.8 (55.9)] 54.9 [77.6 (54.3)] 55.8 [78.5 (55.0)]

2.52 {2.2−5.8}e 9.76 {10.3f, 5.5g} 22.0 {19.6f, 16.4g} 35.2 {31.7f, 29.5g} 45.9 {41.7f, 39.2g} 56.0 {50.8}f 56.5, {51.4}f 55.3 55.0 {48.3}f 56.2 {50.9}f

0.004 [0.006] 0.60 [0.36] 0.18 [0.05] 0.09 [0.04] 0.20 [0.05] 0.11 [0.03] 0.35 [1.30] 0.69 [0.33] 0.53 [0.28] 0.16 [0.11]

0.004 [0.006] 0.54 [0.36] 0.16 [0.04] 0.17 [0.07] 0.20 [0.13] 0.27 [0.13] 0.48 [0.12] 0.76 [0.77] 0.48 [0.09] 0.13 [0.09]

0.004 0.04 0.13 0.11 0.16 0.20 0.44 0.49 0.27 0.13

a The binding energy is defined as the energy sum of all fragments (n × Ewater + Esub) minus the total energy. The PBE-D3 values are in brackets in which the scaled values are in parentheses (scale factor = 0.7). The reference values are in braces. RMSD denotes the deviation in water geometry from the pristine water cluster. bThe value is shown only for the lowest energy orientation. cReferences 11−14. dReferences 11, 12, 14, 19. e References12, 16. fDang−Chang model potential II of ref 9. gPBE/PAW calculation of ref 10.

The ring hexamer adsorbs with little geometry change as long as the hexagon plane in the initial geometry is oriented more parallel to the substrate surface. It converts to the book hexamer when the hexagon plane in the initial geometry is oriented more parallel to the surface normal. (iii) The optimized book hexamer on a graphene-like surface has either open or flipped geometry. The energy difference of the two geometries is almost negligible, and the unfolding angle of the book in the adsorbed state is larger than in the isolated state. The root-mean-square deviation (RMSD)45 of each calculated water cluster geometry with substrates from the corresponding pristine water cluster geometry is reported in Table 2 for (H2O)1−6 on the three graphene-like surfaces by using the Kabsch algorithm46 to superimpose the two geometries. The RMSD values for n = 1 and for the planar ring geometries of n = 3−5, 6 (ring) are less than 0.2 Å. The RMSD value of n = 2 is relatively high as the Hd’s point to different directions from the isolated dimer in the cases of Cor and DBC, but not high in the case of graphene. It arises because the Cor and DBC cases show nonuniform potentials on the molecular surface due to edge H atoms, while the graphene surface shows almost flat uniform potential surface. The RMSD values for nonplanar hexamers (book, bag, cage) are somewhat high due to their facile structure, while that of the rigid nonplanar prism hexamer is small. Binding energies of (H2O)1−6 on the three graphene-like surfaces are also summarized in Table 2. Overall, the binding energies on DBC or graphene are slightly larger than those on Cor, as the H−π interaction increases with the larger substrate. The binding energies of (H2O)2−6 are larger than those reported by Karapetian et al.,9 Leenaerts et al.,10 or Lin et al.17 (Lin et al.’s values are not listed in Table 2 because the substrates are different). This is probably because their calculation methods underestimate the H-bond interactions in water clusters. In the case of water hexamer, the bag, the cage, and the prism structures are not 2D structures, but they are no longer more stable than or nearly isoenergetic to the global minimum energy structure of either ring or book shape. 3.2. Cor/(H2O)1−6/Cor. Interaction energies of (H2O)1,2 as a function of Cor−Cor distance (dz) are plotted in Figure 3. Two kinds of interaction energies are defined as follows. E1 = Etotal − 2 × Esub − n × Ewater

(1)

E2 = E1 − Esub − sub

(2)

Figure 3. Interaction energy as a function of the Cor−-Cor distance (dz) for (H2O)1,2 between a pair of sandwiched Cor molecules, optimized with the SCC-DFTB-3D method. In the (H2O)2 case, the plots are for the dHdd conformer as the initial geometry.

Here Etotal, Esub, Ewater, and Esub-sub are the total electronic energy, the electronic energy of the substrate, the electronic energy of a single water, and electronic energy of the substrate pair at a given layer displacement (dx, dy, dz), respectively. Figure 3 displays the energy curves for a sandwiched substrate pair, i.e., (dx, dy) = (0, 0). Both E1 and E2 curves resemble typical molecular interaction curves; as dz is reduced, the energy decreases slowly to the minimum and then increases rapidly due to the hard wall repulsion. As an H2O monomer is pressured by reducing dz, the E2 minimum is reached when dz = 5.9 Å. Until this point, the geometry of H2O monomer is little changed from that on a single Cor in which the water molecule stands upright while maintaining C2v symmetry (“dd” geometry, Figure 4a). This structure is maintained until dz is shortened to that (5.5 Å) of the E1 minimum. From this point, H2O is optimized with the upper Cor displaced along x or y directions, each time by 0.1 Å and up to ±0.5 Å. In either direction, the energy increases; therefore, the exactly stacked Cor pair provides an optimal environment to stabilize the H2O monomer. As dz is further shortened, the two water hydrogens slowly turn to point horizontally (Figure 4b). We also tried to optimize an initial geometry in which two hydrogens point in opposite directions: one toward bottom Cor, and the other toward top Cor when the two Cor molecules are separated by 9.0 Å. However, this “ud” geometry does not stay stable at the center of the Cor ring. H2O is moved to the Cor edge E

DOI: 10.1021/acs.jpcc.6b06076 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 4. Geometry of H2O monomer between a pair of Cor molecules optimized with the SCC-DFTB-3D method: (a) lowest energy state at (dx, dy, dz) = (0.0, 0.0, 5.5) Å, (b) overcompressed state at (dx, dy, dz) = (0.0, 0.0, 5.0) Å. In the top view of each case, the top Cor is removed for a clearer view.

Figure 5. Geometry of (H2O)2 between a pair of Cor molecules optimized with the SCC-DFTB-3D method: (a) geometry at (dx, dy, dz) = (0.0, 0.0, 7.5) Å, (b) final optimized state at (dx, dy, dz) = (−0.4, −0.3, 5.6) Å. The top Cor and the bottom Cor are shown with blue and black wires, respectively.

when dz is reduced to 8.5 Å. As dz is further reduced, the optimized H2O is placed farther from the edge while the overall energy keeps decreasing. At a fixed dz, the energy of H2O at the edge of the Cor pair is always lower than that of the intercalated H2O. This arises because the energy of Cor dimer plus the energy of H2O is lower than the energy of Cor/H2O/Cor where H2O is intercalated. Namely, the state of intercalation is a local minimum which can be generated from a well-defined precursor state. When the “ud” geometry was optimized inside the DBC pair or the graphene pair, it converted to “dd” geometry. The geometry of (H2O)2 varies in a more complicated fashion as dz is decreased. When the “dHdd” conformer of Cor(H2O)2 is pressured by the top Cor, the adsorbed geometry in Figure 2 is qualitatively maintained until dz reaches 9.0 Å. Then, at shorter dz, the H-bond accepting water is lifted as its dangling OH is pulled by the top Cor. During this process, the whole H2O dimer moves to the edge of the Cor plane (Figure 5a). As dz is reduced to 7.0 Å, the lifted H2O comes back to the same height as the other H2O. In this way, the geometry of the E1 minimum is reached when dz = 5.6 Å. At this stage, the Hd of the H-bond donating water points sideways, and the two Hd’s of the H-bond accepting water point upward and downward. Upon further optimization of the structure by displacing the top Cor along the x and y directions, the two water molecules are translated slightly with a somewhat large change in the pointing directions of the three Hd’s. The final optimized geometry is found at (dx, dy, dz) = (−0.4, −0.3, 5.6) Å (Figure 5b). This geometry will be denoted as “hHud” even though the directions are intermediate between horizontal and up/down. When the “uHdd” conformer of Cor(H2O)2 is pressured by the top Cor, the geometry change undergoes a different route. In this case, the H-bond donating water is lifted as its Hd is pulled by the top Cor when dz is between 7.0 and 9.0 Å. In fact, the energies of these geometries in this range of dz’s are lower by 0.1−1.4 kcal/mol than the corresponding geometries optimized from the “dHdd” conformer. Then, at shorter dz’s, the optimized geometry from the “uHdd” conformer is merged into the geometry from the “dHdd” conformer with the reflection symmetry. At the PBE-

D3/TZVP level, the water dimer tends to stay more toward the central part of Cor. Interaction energy curves for (H2O)3−5 upon changing interlayer distance dz behave similarly to the case of (H2O)1−2. For these sizes of water clusters, the ring structures are known to be stable not only in the isolated phase but also on the graphene surface.10,17 Initially as a planar structure centered on the bottom Cor at dz = 20 Å, all the ring clusters undergo rather a big geometry change when dz is in the range between 7 and 9 Å. In the case of a water trimer, one water molecule whose Hd is pointing up is pulled toward the top Cor at dz = 9 Å and returns to its planar position at dz = 7 Å. In the case of a water tetramer, two water molecules whose Hd’s are pointing up are pulled toward the top Cor at dz = 9 Å and return to their planar position at dz = 7 Å. When dz is reduced to that (6.2 ± 0.5 Å) of the energy minimum, the geometry (in terms of O atoms) is a perfect square. The water tetramer undergoes the geometry change in the following fashion: planar square (>9 Å) → bent square (7−9 Å) → planar square (6−7 Å) → planar rectangle (5−6 Å). A water pentamer behaves similarly. It undergoes the geometry change in the following fashion: planar pentagon (>9 Å) → puckered pentagon (7−9 Å) → planar regular pentagon (6−7 Å) → slightly distorted pentagon (9 Å) → armchair form (7−9 Å) → planar hexagon (