Static Density Functional Study of Graphene–Hexagonal Bilayer Ice

Mar 18, 2014 - Previous modeling relevant to the graphene–ice-1h interaction was ... and librations, and cellular ΔH – TΔS offers at least a bal...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/JPCA

Static Density Functional Study of Graphene−Hexagonal Bilayer Ice Interaction David J. Anick Laboratory for Water and Surface Studies Department of Chemistry, Tufts University, 62 Pearson Avenue, Medford, Massachusetts 02155, Unites States S Supporting Information *

ABSTRACT: Periodic static ab initio studies are conducted of hexagonal bilayer ice (HBL) and basal layers of ice-1h adsorbed on graphene using the model BLYP-D in CRYSTAL09. Eight highsymmetry periodic forms of HBL are optimized, of which four have lower energy; their electronic binding energy to graphene is ∼1.6 kcal/mol per abutting H2O. Optimized geometries have the property of maximizing the occurrence of a certain O−H−C alignment motif. One lattice is selected for more detailed study. Its 2-D shear translation potential energy surface is found to have barrier heights in two zigzag directions of ∼140 cal/mol per abutting H2O. A second hexagonal bilayer can be added and the electronic binding energy drops from ∼1.7 to ∼1.0 kcal/mol per abutting H2O. For ice-1h monolayer adsorbed on graphene, a proton-ordered form in which half of the O’s nearest the graphene carry a proton pointing toward graphene is preferred over protonordered forms in which either all or none of those O’s have H’s pointing toward graphene. Cohesive energy for two-layer ice-1h on graphene is 0.66 kcal/mol of H2O higher than for HBL, supporting experimental evidence that the graphene+HBL isomer is more stable. However, the HBL and two-HBL structures are unstable or at best metastable with respect to four layers of ice-1h. hexagonal tessellations (its “faces”) in register, with O’s at the vertices, and one H-bond along each hexagon edge (called “facial” H-bonds) and one H-bond linking each O to its counterpart O in the opposing plane (called “joining” Hbonds). Thus each O is 4-coordinated but its local geometry is distorted or buckled from tetrahedral, its nearest neighbors sitting at four of the five positions of a tripyramidal arrangement (cf. Figure 1a). Waters that serve as donor in a joining bond are “joining” H2O’s; otherwise they are “facial”. The study reported herein used static 2-D periodic ab initio modeling to explore Kimmel’s HBL-on-graphene structure and to compare it with ice-1h on graphene. The Pt(111) backing has a “negligible” effect on graphene’s properties, so modeling it as free-standing graphene is reasonable.21 The HBL structure has an interesting history, starting with its first discovery by Koga18 as a crystalline phase in a confined geometry. Molinero et al.16 conducted dynamical modeling of HBL and several other confined phases having two to four molecular layers and found that HBL, and only HBL, has an anomalously high melting temperature well above 0 °C. Jinesh22 found evidence for ice at room temperature on graphene during atomic force microscopy measurements and posited HBL formation as the most likely explanation. Kirov23,24 computed the residual entropy of HBL, finding it to be just 43% of that of a two-layer ice-1h. In ref 25 it was

1. INTRODUCTION Since the development of an efficient exfolation method for producing high quality graphene sheets by the Manchester group in 2004,1,2 the study of graphene’s remarkable electrical, thermal, and mechanical properties3,4 has grown rapidly. Interactions of graphene with water and ice are important to understand in view of the ubiquity of humidity, graphene’s potential as a coating for electrodes, e.g., in solar cells,5 and its promise as a material for water treatment.6 As a polyaromatic without polar side groups, graphene is normally considered hydrophobic, but Feller and Jordan demonstrated as early as 2000 that a single H2O would orient with one O−H bond toward the sheet, and its binding to graphene would be stronger than the binding in the water dimer.7 Bermudez and Robinson8 subsequently obtained the value of 207 meV (4.77 kcal/mol) for single H2O binding. Further evidence that graphene is no ordinary hydrophobe comes from observations that the contact angle between water and graphene varies greatly with conditions, and is tunable;9 and that desorption follows a complex 3-peak pattern reflecting competition between water−water bonds and water-graphene bonds.10 Kimmel et al.11 discovered that under certain conditions, a thin water film growing on graphene on Pt(111) adopts a hexagonal bilayer structure with four H-bonds at each O rather than a more conventional tetrahedrally coordinated ice-like structure. The context for this surprising result was a series of experimental and modeling studies of nanoconfined water, some of them using graphene for the nanopore walls, in which exotic water phases were reported.12−20 The hexagonal bilayer (henceforth HBL) structure consists of two parallel planar © XXXX American Chemical Society

Special Issue: Kenneth D. Jordan Festschrift Received: January 13, 2014 Revised: March 17, 2014

A

dx.doi.org/10.1021/jp500360n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Figure 1. (a) Hexagonal bilayer (HBL); (b) HBL on graphene; (c) ice-1h monolayer on graphene. Gray = C, white = H. For (a), blue = joining O’s and red = facial O’s. For (b) and (c) red = O. Selected BLYP-D-optimized H-bond lengths and O−O−O angles are included. Orange dummy atoms indicate the corners of a rectangular crystallographic cell for (c). Orange and yellow dummy atoms together define the primitive cell vectors and give the corners of a rhomboidal primitive cell for (a).

graphene−water and found values ranging from 88 to 138 meV per H2O for the binding energy. Closest approach of graphene to any water O ranged from 3.09 to 3.38 Å for the vdW-corrected functionals. Questions that this study considered included the following. • What is the binding energy of HBL on graphene, and what is the alignment in the optimal geometry? Do Hbond directions in the HBL make much of a difference? • What are the barrier heights for shear motion of the HBL along a graphene sheet? • What is the binding energy for 1 or more basal ice-1h layers on graphene and is there a preference for the basal ice surface to abut the graphene with non-H-bonded H’s, with oxygen lone pairs, or with some of each? • What are the energetics of HBL vs 2L1h (i.e., a double basal layer of ice-1h), either free or adhered to graphene? What temperature effect would be predicted? • What can be said about the relative favorability of 2HBL (two stacked HBL bilayers) vs 4L1h (four layers of ice1h), either free or adhered to graphene?

noted that hexagonal bilayers can be stacked to form a 3-D phase, “hexagonal bilayer water” (HBW). In HBW, H-bond directions have a much greater impact on structural parameters than for tetrahedral ices. It was noted that eight distinct protonordered isomers of HBL can be described that have high symmetry; i.e., the entire lattice is generated by just two H2O’s. Optimal stacking arrangements were most often monoclinic lattices and appeared to be guided by the preference for a motif in which an O−H---O makes an angle of 108° ± 10° whereas the O−O separation is between 2.95 and 3.35 Å. Due to many (perhaps coincidental) similarities between HBW and high density amorphous ice (HDA), the question was raised of whether islands of HBW might be present to a significant extent in HDA. Graphene has the remarkable property that its bond length of 1.42 Å26,27 allows it to approximate closely the dimensions of both free-standing HBL and a basal ice-1h layer. The alignments of graphene with HBL and with a monolayer of ice-1h are shown in Figure 1b,c. Though ice normally has Hbonds of length 2.76 Å, when templated by the graphene as in Figure 1b, the facial H-bonds become 2 × 1.42 = 2.84 Å; i.e., they are stretched by about 3%. Because the basal layer of ice1h is puckered, its projection onto the basal plane consists of hexagons whose side length is ideally (8/9)1/2 × 2.76 = 2.60 Å. When templated by graphene, this shrinks to √3 × 1.42 = 2.46 Å, a decrease of ∼5%. The decrease is achieved by a combination of shortening the bond and increasing the puckering of the basal surface. Previous modeling relevant to the graphene−ice-1h interaction was undertaken by Leenaerts et al.,28,29 who explored binding of (H2O)n clusters, 1 ≤ n ≤ 6, to graphene. Their largest clusters, (H2O)6, consisted of a ring hexamer, whose optimized adherent geometry closely resembles a hexagon of Figure 1c. The tendency for a single H2O on graphene to adopt an orientation resembling H2O donating in an H-bond initially made some think this would be water’s preferred setup in larger groupings also, but Leenarts found that larger clusters had both up- and down-pointing free H’s. Politano’s spectroscopic work with ice monolayers on graphene likewise suggested some free H’s are not interacting with the graphene;21 hence, they presumably are pointing away from it. In a recent work Ambrosetti et al.30 obtained binding energies for a single or double basal ice-1h layer on graphene for several dispersioncorrected and uncorrected functionals. They noted how the dispersion correction is essential in any DFT study of

2. METHODS All reported calculations were done using the two-dimensional SLAB option in CRYSTAL0931 running on a PQS Quantum Cube.32 The vacuum between slabs is 500 Å. Density functional theory (DFT) calculations were done using the Becke−Lee− Yang−Parr functional33 with the Grimme dispersion correction.34,35 Many-body interaction terms would improve the accuracy of the dispersion correction36−40 but only DFT-D is implemented in CRYSTAL09. On H and O the triple-ζ Ahlrichs basis set (TZVP) was used, and on C, the double-ζ set (DZVP).41 Large basis sets can cause problems with linear dependence for periodic systems with atoms closely spaced, and this was true in trials of TZVP for graphene in CRYSTAL09. BLYP-D has undergone extensive benchmarking in studies of water,42−44 and although it does a good job of reproducing many of water’s properties, it is also known to overbind ice-Ih to some extent. All reported optimizations included optimization of cell parameters and included computation of frequencies. CRYSTAL09 only computes “cellular” vibrational modes; i.e., a Hessian is computed numerically for 3N degrees of freedom if there are N atoms in a primitive cell. The harmonic approximation is then used for frequencies. CRYSTAL09’s B

dx.doi.org/10.1021/jp500360n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Figure 2. Hexagonal prism fragments of the eight high-symmetry lattices, showing BLYP-D-optimized geometries of the lattices when adherent to graphene. Gray = C, white = H, red = O of facial H2O, blue = O of joining H2O. To generate the full lattice, translate in the X direction so that the left edge of the copy coincides with the right edge shown, and translate in the Y direction so that the lowermost O of the copy is in an H-bond with the uppermost O shown. “Z” means that the joining O’s (or the facing O’s) form zigzag lines, whereas “T” means that the joining O’s (or the facing O’s) form a triangular sublattice. “A” means the direction of every H-bond is antiparallel to that of the bond above or below it, whereas “B” means that both parallel and antiparallel can occur. “P” means that within each hexagon, opposite sides have their H-bond directions parallel, whereas “M” means that for opposite side bonds both parallel and antiparallel relationships occur (“mixed”). The translation to the scheme of ref 25 would be T↔F, Z↔J; A↔W, B↔E; P↔E, M↔W.

Table 1. BLYP-D Cell Parameters and Energies for Free-Standing Hexagonal Bilayer Isomers and Two-Bilayer Isomers (Cohesive Energies in kcal/mol) .

a

cell parameters

cohesive energies 0

geom

NC

NW

sym group

a

b

E

E0 + cZPE

cΔG298

ZBP ZBM ZAP ZAM TBP TBM TAP TAM 2TAPa 2TAPb 2TAPc

0 0 0 0 0 0 0 0 0 0 0

8 8 8 8 8 8 8 8 16 16 16

P2an P21/b11 P21/b11 P2an C211a P21/b11 C2/m11a P2/b11 C2/m11a C2/m11a Cm11a

4.761 8.140 8.134 4.791 4.773 4.759 4.638 4.800 4.599 4.660 4.639

8.131 4.752 4.761 8.073 8.090 8.109 8.328 8.041 8.371 8.314 8.340

−135.23 −135.37 −135.58 −135.27 −135.82 −135.99 −136.22 −135.87 −277.87 −278.28 −278.42

−109.62 −110.00 −110.01 −109.72 −110.26 −110.54 −110.69 −110.15 −226.69 −226.92 −227.10

−25.54 −25.94 −25.94 −25.65 −26.22 −26.43 −26.56 −26.07 −62.15 −61.79 −62.39

Geometries with symmetry C211, C2/m11, or Cm11 were computed using supercells with symmetry P211, P2/m11, or Pm11, respectively.

cellular figures would be more likely to cancel out. This was achieved by using supercells where necessary, for the lattices with greater symmetry. Tables report E0 (electronic energy), E0 + cZPE, and cellular free energy at T = 298 K, denoted cΔG298. Studies indicating that these thin layers can remain crystalline at room temperature16,22,45,46 make 298 K a reasonable choice of reference temperature for this work. How much difference supercells make, was benchmarked by computing one structure using its primitive cell as well as with supercells consisting of 2, 4, 6, or 8 primitive cells: see the Supporting Information.47 For a periodic structure with NC carbons and NW waters in the unit cell, its “cohesive energy” is obtained by subtracting

zero-point energy is based on these 3N modes, so long-range intercellular phonon modes are omitted, and for clarity we will use the notation cZPE. Cellular modes were likewise applied for enthalpy and entropy calculations, which were used in finite temperature results. How valid are cZPE and cellular entropy and enthalpy is a fair question, but the calculation does verify a local minimum or a transition state. It provides approximate vibrational frequencies for O−H stretches, H−O−H bends, and librations, and cellular ΔH − TΔS offers at least a ballpark for the finite temperature correction term. Wherever possible, comparisons only were made among lattices that had the same number of H2O’s per cell, so that any errors deriving from the C

dx.doi.org/10.1021/jp500360n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Table 2. BLYP-D Cell Parameters and Energies for Graphene-Adherent Hexagonal Bilayer Isomers and Two-Bilayer Isomersa cell params geom + gra (Gra) ZBP ZBM ZAP ZAM TBP TBM TAP TAM 2TAPa 2TAPb 2TAPc 2TAPd

NC 16 16 16 16 16 16 16 16 16 16 16 16 16

NW

sym gp

0 8 8 8 8 8 8 8 8 16 16 16 16

b

P6mm Pb11 Pb11 Pb11 Pb11 P1c Pb11 Cm11b Pb11 Cm11b Cm11b Cm11b Cm11b

cohesive energies

binding energies

a

b

E0

E0 + cZPE

cΔG

E0

E0 + cZPE

cΔG

4.971 8.571 8.571 8.571 8.570 4.956 4.951 4.954 4.952 4.938 4.941 4.940 4.940

8.613 4.959 4.959 4.959 4.959 8.576 8.585 8.579 8.584 8.563 8.561 8.565 8.565

0.00 −141.40 −141.38 −141.79 −141.41 −142.21 −142.36 −143.20 −142.30 −281.12 −282.57 −282.47 −282.49

0.00 −116.25 −116.29 −116.68 −116.34 −117.55 −117.38 −117.36 −117.51 −230.38 −231.67 −231.73 −231.89

0.00 −35.75 −36.02 −36.20 −36.00 −37.16 −38.48 −38.41 −38.59 −72.68 −72.72 −73.17 −73.59

1.54 1.50 1.55 1.53 1.60 1.59 1.74 1.61 0.81 1.07 1.01 1.02

1.66 1.57 1.67 1.65 1.82 1.71 1.67 1.84 0.92 1.19 1.16 1.20

2.55 2.52 2.56 2.59 2.73 3.01 2.96 3.13 2.63 2.73 2.69 2.80

The first row refers to graphene with no water. Cohesive energies are in kcal/mol. and binding energies are in kcal/mol of H2O(a). bComputed as Pm11. cOblique lattice, γ = 89.97°.

a

from its energy, the quantity NCEC + NWEW, where EC is the per-carbon energy of graphene and EW is the energy of gas phase monomeric H2O. For each lattice we report three binding and cohesive energies, i.e., for E0, E0 + cZPE, and cΔG298. [The three values of EC for this model are (in au) −38.105734, −38.099611, and −38.099682. The corresponding values for monomeric H2O are −76.443145, −76.422758, and −76.399605.] Units of “kcal/mol” mean per mol of unit cells unless “mol of H2O” is specified in which case it has been divided by NW. “Binding energy” refers to the absolute value of the interaction energy of a unit such as a monolayer or bilayer with graphene. When presented as “per mol of H2O(a)”, cellular energy has been divided by NA, the number of “abutting” H2O’s, i.e., water molecules of the unit cell in the single layer closest to graphene. Because the dispersion force varies as O(r−6), we reason that abutting H2O’s are responsible for nearly all of the graphene−water interaction, so the average over NA may be better for comparisons than the average over NW. Binding energy provides insight about a particular structure whereas cohesive energy is appropriate for thermodynamic comparisons of different structures.

other HBL isomers, its cell dimensions are notably elongated in the armchair direction. The cell aspect ratio for TAP is 1.80, contrasted with 1.68−1.71 for the other seven. Perfectly regular hexagons would produce a cell aspect ratio of √3 ≈ 1.732. TAP is unique in having all its facial H2O’s align their dipole moments either parallel or antiparallel to the same armchair axis (vertical in Figure 1a,b). In planar hexagons, the angle between two facial H-bonds that share a donor O is strained (120° vs the gas phase monomer angle of 104° or the tetrahedral angle 109.47°). Free-standing TAP reduces this angle to 116.2°, which TAP can do because all facial H2O’s undergo the same strain reduction. In the other HBL lattices, reducing one H2O’s O−O−O angle would come at the cost of increasing the O−O−O angle for some other H2O’s. All of the “T” isomers have lower energy than all of the “Z” isomers. The mean difference for E0 (E0 + cZPE, cΔG298) is 0.62 kcal/mol (0.57, 0.55). In the T form, nearest-neighbor joining bonds are antiparallel, whereas in Z isomers some joining H2O’s bond to other joining H2O’s and those adjacent joining bonds are parallel. When graphene and HBL are optimized together, there is much less variation in the cell dimensions (cf. Table 2). Cell aspect ratios range from 1.728 to 1.734. Mostly, the water accommodates to the graphene, but to a slight extent graphene’s dimensions also adjust to the presence of the water. Cohesive energy is consistently lower for the graphenebound T isomers than for the Z isomers. The T-vs-Z difference is ∼1 kcal/mol for E0 and E0 + cZPE and ∼2 kcal/mol for cΔG298. The strength of the interaction per abutting H2O is about one-third that of an H-bond. Figure 2 compares how the eight lattices align on the graphene in their lowest energy geometries. Looking first at the T isomers, the O’s of facial H2O’s sit (very close to) just above the center of a graphene hexagon. The O’s of joining H2O’s sit just above a C with the O−H that is part of a facial bond lying over part of a C−C bond. The Z isomers are not able to distribute their facial and joining and H2O’s in this way, and it is a reasonable speculation that their higher energies result in part from being frustrated from achieving the more ideal setup that is available to the T isomers. There may be a more fundamental consideration behind the optimized geometries. In ref 25 it was noted that minimum energy geometries for adjacent bilayers in HBW typically

3. RESULTS AND DISCUSSION A. Proton Ordered Hexagonal Bilayers. Eight “high symmetry” proton ordered isomers of HBL were studied. They are classified by a set of three binary letter choices, (T or Z)(A or B)(P or M), e.g., the HBL of Figure 1a is TAP. Figure 2 shows an hexagonal prism fragment of each lattice. (See caption of Figure 2 for an explanation of the notation.) They are the building blocks for the fourteen orthorhombic and monoclinic forms of HBW.25 Each was optimized both independently and adherent to graphene. On graphene the unit cells consisted of C16(H2O)8. Table 1 provides data for the free-standing HBL’s and Table 2 for HBL’s adsorbed onto graphene. The tables also give the crystallographic cell dimensions: the crystallographic a and b axes can correspond to the zigzag and armchair axes of the graphene, or vice versa. Standardized representations of the 80 layer groups in crystallography48 determine which axis is a and which is b. Crystal coordinates for all structures listed in Tables 1 and 2 are provided in the Supporting Information.47 Among the free-standing bilayers, lowest energy occurs for the TAP isomer (Figure 1a), which is also the isomer that generates the lowest-energy forms of HBW.25 Compared to the D

dx.doi.org/10.1021/jp500360n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

preventing it from coming as close to the O−H---O as O’s come in HBW. B. Shear Potential Energy Surface (PES) and TwoBilayer Study. Shearing of an ice layer along graphene, i.e., translation parallel to the sheet, was considered by Ambrosetti et al.30 They found the barrier to sliding an ice monolayer was 2 orders of magnitude smaller than for graphene−graphene shear (410 meV). They also obtained the low barrier of 42 meV for shear in the presence of water intercalated between graphene layers, suggesting a mechanism for graphene’s well-known use as a lubricant when some humidity is present. In ab initio and classical MD studies of bulk water adjacent to graphene, Rana and Chandra49 discovered that the immediately abutting monolayer displayed slowed dynamics suggestive of partial freezing while also permitting faster-than-bulk translational relaxation along the sheet, consistent with a low shear barrier. Reference 25 discussed a low-barrier “sliding PES” for stacked hexagonal bilayers sliding relative to one another. Six of the eight HBL isomers have two local minima for how they can stack, and the low barrier suggested that 3-D HBW would occur as a quasi-crystal that combines both offsets in a stochastic mix. We computed both the barriers to HBL sliding along graphene and the energetics of adding a second HBL to Gra +HBL. The TAP isomer was used for both studies. A 2-D vertical PES was computed starting from the optimized geometry of TAP on graphene. The water was translated relative to the graphene, in steps of 0.1 fractional coordinates, along the two primitive cell vectors. The resulting plot is given as Figure 4. The minima are at (0, 0), (0.5, 0), (1, 0), (0, 0.5), etc. because periodicity of the graphene at twice the spatial frequency of the lattice means that the PES inherits the smaller periodicity. Barriers in the primitive cell vector directions, which are zigzag directions for both the graphene and the HBL, are ∼140

exhibit a motif in which a facial O−H of one bilayer aligns with an O of the other bilayer in such a way that the O−H---O angle is 108° ± 10° and the O−O distance is 315 ± 20 pm. Likewise, in every optimized graphene−HBL geometry we find that each facial O−H adjacent to graphene participates in one or more setups we call the “OHC motif”. The OHC motif occurs if the O−H---C angle is 105° ± 10° and the O−C distance is 328 ± 12 pm (cf. Figure 3). This is not true of randomly arranged

Figure 3. Fragment of the right-hand edge of ZBM in Figure 2, illustrating the relationships in the OHC motif. O1−H1−C1 = 108.0°, O1C1 = 3.27 Å. O2−H2−C2 = 99.8°, O2C2 = 3.19 Å. O2−H3− C3 = 106.3°, O2C3 = 3.38 Å.

bilayers nor of intermediates in the optimization process. A possible hypothesis is that anisotropy of the electron density around an O−H---O hydrogen bond permits close approach by a fourth atom (O or C) only when that fourth atom approaches from a characteristic direction relative to the O−H. The vdW interaction is then maximized when the fourth atom enters this low electron density pocket. We emphasize that this would only apply to O−H if the H is already in an H-bond with another O. Behavior of free O−H’s would be dominated by electrostatics and they would not have a tendency to form an OHC motif. The overall size of the OHC effect is small, perhaps in part because each C is constrained by its strong bonds to other C’s,

Figure 4. Vertical potential energy surface obtained by shear translation of water relative to graphene, for the Gra+TAP system. Horizontal units are fractional coordinates in the primitive cell (cf. Figure 1a), indicating relative motion along two zigzag axes. Vertical units are cal/mol of H2O(a). E

dx.doi.org/10.1021/jp500360n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Table 3. BLYP-D Cell Parameters and Energies for Free and Graphene-Adherent Ice-1h Layersa cell parameters

cohesive or (binding) energies

geometry

NW

NC

sym gp

a

b

γ

1L1h(m) 2L1h 3L1h 4L1h Gra+1L1h(T) Gra+1L1h(A) Gra+1L1h(m) Gra+2L1h Gra+3L1h Gra+4L1h

4 8 12 16 4 4 4 8 12 16

0 0 0 0 12 12 12 12 12 12

P1 P1 P1 P1 Cm11b Cm11b P1 P1 P1 P1

4.483 4.455 4.437 4.429 4.312 4.309 4.307 4.310 4.313 4.315

7.647 7.637 7.635 7.636 7.471 7.468 7.460 7.464 7.468 7.473

89.62 89.95 90.07 90.12 90.00 90.00 89.99 90.00 90.01 90.03

E0 −53.63 −129.12 −205.47 −281.94 −58.21 −55.97 −62.19 −137.94 −214.11 −290.34

E0 + cZPE

(2.14) (2.20) (2.16) (2.10)

−43.35 −105.97 −169.61 −233.25 −47.93 −46.02 −51.86 −114.61 −178.66 −241.46

(2.13) (2.16) (2.26) (2.05)

cΔG298 −2.33 −23.66 −46.60 −69.83 −8.66 −6.67 −12.50 −34.65 −56.97 −77.19

(2.54) (2.75) (2.59) (1.84)

a

Cohesive energies in kcal/mol, binding energies (in parentheses) in kcal/mol of H2O(a). (T) = all-toward isomer, (A) = all-away, (m) = 50-50 mixed. bComputed as Pm11.

where the number of adjacent pairs of daa waters would be minimized in the monolayer. The last criterion is known to be important when low energy isomers of 3-coordinated water clusters are sought.50,51 Hirsch and Ojamäe classified 16 proton-ordered isomers of ice-Ih.52 The only isomer in their list that fits all of these specifications is their isomer #14. The 50− 50 monolayer and all multilayer forms of ice-Ih studied herein are slabs of the #14 isomer. It is illustrated on graphene in Figure 1c. We did not attempt a comprehensive study of the effects of proton ordering as we did for the HBLs. For depictions of the “all toward” monolayer on graphene (Gra +1L1h(T)) and the “all away” version (Gra+1L1h(A)); please see Figure 6 in the Supporting Information.47 We were unable to optimize the free ice-XI monolayer. It has a strong dipole moment and tends to crumple during optimization, ending up as something rather different from a flat water monolayer. The 50−50 monolayer does have a local minimum. Energies of ice layers with and without graphene are compared in Table 3. The binding energy for the 50−50 monolayer is 2.14 kcal/mol of H2O(a) for E0 (2.13 and 2.54 for E0 + cZPE and cΔG298, respectively) and 2.20 (2.16 and 2.75, respectively) for the ice-1h bilayer. These figures are consistent with those of Ambrosetti et al., which ranged from 2.03 to 2.80 kcal/mol for the vdW-corrected functionals.30 Comparing cohesive energy absolute values for ice-1h monolayers on graphene, we see that 50−50 > all-toward > all-away. Viewing cohesive energy a function of the percentage of free H’s that point toward graphene, it can be approximated as a quadratic curve passing through (for E0) the points (0, −55.97), (50, −62.19), and (100, −58.21). For E0, as well as for E0 + cZPE and for cΔG298, the minimum occurs at 55% toward. This is consistent with the emerging picture that both toward and away orientations occur but there is a slight preference for toward. It might be possible to determine this percentage experimentally by saturating the away-pointing H’s of Gra +1L1h with trimethylamine (TMA), and then measuring the ratio of TMA to water recovered.53 As seen in the last four rows of Table 3, binding energy decreases with increasing number of ice layers, except for the anomaly of an increase on going from one to two layers. The effect is much less dramatic than for the addition of a second hexagonal bilayer on Gra+HBL. We suspect this is because HBL has flat faces and must accommodate by altering its Hbond lengths, whereas ice-1h can accommodate by altering some O−O−O angles, which is less energetically costly. This explanation is supported by the small but steady increase of

cal/mol of H2O(a). This means there is a pathway that follows a slightly wavy course but arrives at the next periodic minimum after cresting a saddle 140 cal/mol of H2O(a) above Figure 1b. This figure is comparable to Ambrosetti’s (∼100 cal/mol) for the ice monolayer barrier. The “hills” of Figure 4 peak at 390 cal/mol of H2O(a) above the minimum. Integrating the whole PES over the unit cell gives an average height of 180 cal/mol of H2O(a). We identified four local minima among ways a second copy of the TAP bilayer can be stacked atop Gra+TAP. For two of them, the orientation of the upper bilayer is the same as the lower. Looking for orientation-preserving ways of stacking HBLs to make HBW, ref 25 found two inequivalent monoclinic angles that yielded local minima for stacked TAPs, denoted the “a” and “b” forms, so we call the corresponding 2TAP isomers “2TAPa” and “2TAPb”. Two other local minima reverse the upper copy of TAP in the b direction relative to the lower bilayer, and the notations are “2TAPc”’ and “2TAPd”. In the absence of graphene the “c” and “d” forms are equivalent, but the symmetry operation that equates them is lost when graphene is on one surface. Figure 5 (Supporting Information47) illustrates the four alignments. Their energetics are included in Tables 1 and 2. The 2TAPa isomer has the smallest binding and there is little difference among the other three isomers. All show a significant increase of binding energy with temperature. The 2TAP isomers have lower binding energy than TAP, which can be understood in terms of bond length strain. When a single HBL accommodates to graphene, the bonds must stretch, and the binding energy represents the net benefit of the graphene interaction minus the bond strain. When two HBLs accommodate as a stacked unit, the binding energy now represents the net benefit of the same graphene interaction minus twice as much bond strain. This can be seen in comparing the cell parameters for Gra+2TAP vs Gra +TAP: the b dimension contracts by about 2% in going from Gra+TAP to Gra+2TAP. Some of the strain energy is held by the graphene, even though, as with single bilayers, our calculations support the general understanding that water does most of the accommodating. C. Ice-1h on Graphene. We did computations for each of three ice-1h monolayers on graphene: one with all free H’s pointing toward graphene, one with all free H’s pointing away, and one with a 50−50 mix. For the first two, a basal layer of iceXI could be used. For the third one, we wanted a protonordered form of ice-Ih that is antiferroelectric, where each basal layer donates equally often in the +c as in the −c direction, and F

dx.doi.org/10.1021/jp500360n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

both the a and the b cell dimensions as the number of layers grows, indicating that graphene shares proportionately more of the strain each time another layer is added. D. Comparisons of Free and Graphene-Adherent HBL and 2HBL with Ice-1h Layers. Cohesive energy of TAP, which is the lowest-energy free HBL, lies below that of the free ice-1h bilayer, denoted 2L1h. Cohesive energy differences “per H2O”, i.e., cellular energies divided by NW, are listed in Table 4.

any Gra+2HBL, there are preliminary reports that experimental systems with two or more HBL bilayers on graphene have been made,55 so the comparison bears a closer look. Returning to Table 2 and binding data, a second TAP bilayer weakens the attachment of TAP to graphene by ∼0.6 kcal/mol of H2O(a), an effect that we have attributed to the cost of stretching the Hbonds of the second bilayer. We have also seen that shearing a bilayer costs a maximum of 0.39 kcal/mol of H2O(a) at a hilltop. The average over the whole shear PES is 0.18 kcal/mol of H2O(a). If an HBL were adjacent to graphene but the positions of the O’s relative to the C’s were uncorrelated in the a and b directions, binding would go down by only 0.18 kcal/ mol of H2O(a). This means the cost of falling out of optimal alignment with the graphene is significantly less than the cost of keeping the lattice’s dimensions. Once a second bilayer is added, Gra+2HBL can lower its total energy if the 2HBL and graphene do not try to operate as a single lattice. If the 2HBL stays adjacent to graphene but adopts its own optimal cell parameters, the abutting O and C positions will become essentially uncorrelated. That will cost (on average) 0.18 in misalignment but will regain the lost 0.6, a net advantage of ∼0.42 kcal/mol of H2O(a). The picture painted here is that graphene epitaxy controls the deposition of an initial single HBL, but once the HBL begins to template a second HBL, it is advantageous for the emerging two-bilayer system to unhook from its mooring and for parts of it to slide along graphene as needed to achieve its own independent optimal geometry. A “two-lattice” geometry cannot be modeled as a periodic system, unless one uses a very large cell that is a least common multiple of both lattices’ dimensions. By our crude calculations, the energy gained would not be enough to overtake 4L1h, but combined with a favorable temperature effect, and taking into account the tendency for BLYP-D to overbind ice-1h, this alternate form of Gra+2HBL might be a stable or metastable configuration under some conditions. Our findings do support the idea that nHBL has a reasonable chance to exist in confined geometries. 2HBL is narrower than 4L1h so a slit nanopore whose width is ideal for 2HBL would force 4L1h to be compressed and distorted, giving the advantage to 2HBL. When Molinero et al.16 explored systems of up to four H2O layers in a slit nanopore with hydrophobic walls, they did not include 2HBL as one of their geometries. Given the high stability and high melting point of the single HBL in a confined geometry, confined systems consisting of nHBL for small n would certainly be interesting to consider.

Table 4. BLYP-D Cohesive Energy Differences (kcal/mol of H2O), for Two to Four layers of stacked TAP vs the Same Number of Layers of Ice-1h, i.e., (E(n/2 TAP) − E(nL1h))/ NW, Either Free or Adherent to Graphene state

n

E0

E0 + cZPE

cΔG298

free free free Gra Gra Gra

2 3 4 2 3 4

−0.89 −0.15 0.22 −0.66 0.11 0.49

−0.59 0.06 0.38 −0.34 0.34 0.60

−0.36 0.18 0.46 −0.47 0.08 0.23

Free TAP beats free 2L1h by 0.89 kcal/mol of H2O for E0 (respectively 0.59, 0.36 for E0 + cZPE, cΔG298). Bond angle strain in HBL is more than compensated by having the full complement of four H-bonds at each O. (Recall that 2L1h has an equal number of 3- and 4-coordinated O’s.) The temperature effect favors 2L1h. For thicker layers, n-layer ice1h (notation nL1h) can be compared against n/2 stacked HBL’s if n is even; or for n odd, against a stacked HBL setup consisting of (n − 1)/2 bilayers with half of its covered area having an additional bilayer (notation is n/2 HBL). For n = 3 the comparison is close and for n = 4 it favors 4L1h (cf. Table 4). Binding to graphene is slightly stronger for 2L1h than for HBL, but not enough to cancel out the intrinsically lower energy of HBL. Gra+TAP is favored over Gra+2L1h by 0.66 kcal/mol of H2O for E0 (respectively 0.34, 0.47 for E0 + cZPE, cΔG298). On graphene the temperature effect favors Gra+HBL. Again the comparison is within the margin of error for n = 3 and favors Gra+4L1h for n = 4. The favorability of Gra+HBL over Gra+2L1h is in line with Kimmel’s findings that Gra+HBL occurs on graphene− Pt(111), with “no confinement necessary”. However, Kimmel et al. also reported that once warmed above 130 K, the system reverted to ice-1h, even though our computed temperature effect would not have predicted that. The calculation may be correct if the experimental system is not converting entirely to 2L1h but instead is generating regions of nL1h for n ≥ 3. Conversion of Gra+HBL to Gra+2L1h shrinks the region of graphene that is covered by 25%, leaving the remainder bare of water. Unless exact coverage is known before and after conversion, the system might leave gaps of more than 25% and bunch up the ice-1h in some places, yielding patches of 3 or more layers. Above 130 K, ice has enough thermal energy to drive some rearrangements that break and reorder H-bonds; e.g., the conversion of HDA to low density amorphous ice occurs around 135 K.54 The conversion of Gra+HBL to Gra +nL1h will be exothermic if the average n for the occupied regions exceeds a cutoff that is between 2 and 4 and is most likely close to 3. E. Speculations on 2HBL. Although our results clearly favor 4L1h over any 2HBL and clearly favor Gra+4L1h over

4. CONCLUSION We have presented here the beginnings of a static density functional modeling of the graphene−hexagonal bilayer ice system. Although we have made reasonable choices for models and methods, we must caution against taking any calculations herein as definitive. BLYP-D is just one vdW-corrected functional and we have seen that there is wide variation among such functionals, in their prediction of water−graphene interaction energies. In addition to dispersion corrections that include many-body terms, inclusion of self-interaction error effects may generally improve the description of noncovalently bonded systems.40 Some uncertainty may also arise when finite temperature effects are inferred from relatively small cells. Still, several take-home messages from our work are worth noting. We found a moderately strong interaction between graphene and HBL. This adds to the ways in which graphene is G

dx.doi.org/10.1021/jp500360n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

(9) Zhou, H.; Ganesh, P.; Presser, V.; et al. Understanding Controls on Interfacial Wetting at Epitaxial Graphene: Experiment and Theory. Phys. Rev. B 2012, 85, 035406. (10) Souda, R. Nanoconfinement Effects of Water on Hydrophilic and Hydrophobic Substrates at Cryogenic Temperatures. J. Phys. Chem. C 2012, 116, 20895−20901. (11) Kimmel, G. A.; Matthiesen, J.; Baer, M.; et al. No Confinement needed: Observation of a Metastable Hydrophobic Wetting TwoLayer Ice on Graphene. J. Am. Chem. Soc. 2009, 131, 12838−12844. (12) Bai, J.; Zeng, X. C. Polymorphism and Polyamorphism in Bilayer Water Confined to Slit Nanopore Under High Pressure. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 21240−21245. (13) Chialvo, A. A.; Vlcek, L.; Cummings, P. T. Surface Corrugation Effects on the Water − Graphene Interfacial and Confinement Behavior. J. Phys. Chem. C 2013, 117, 21875−21886. (14) Giovambattista, N.; Rossky, P. J.; Debenedetti, P. G. Computational Studies of Pressure, Temperature, and Surface Effects on the Structure and Thermodynamics of Confined Water. Annu. Rev. Phys. Chem. 2012, 63, 179−200. (15) Johnston, J. C.; Kastelowitz, N.; Molinero, V. Liquid to Quasicrystal Transition in Bilayer Water. J. Chem. Phys. 2010, 133, 154516 (2010). (16) Kastelowitz, N.; Johnston, J. C.; Molinero, V. The Anomalously High Melting Temperature of Bilayer Ice. J. Chem. Phys. 2010, 132, 124511. (17) Kim, J.-S.; Choi, J. S.; Lee, M. J.; et al. Between Scylla and Charybdis: Hydrophobic Graphene-Guided Water Diffusion on Hydrophilic Substrates. Sci. Rep. 2013, 3, 2309. (18) Koga, K.; Zeng, X.; Tanaka, H. Freezing of Confined Water: A Bilayer Ice Phase in Hydrophobic Nanopores. Phys. Rev. Lett. 1997, 79, 5262−5265. (19) Koga, K.; Tanaka, H. Phase Diagram of Water Between Hydrophobic Surfaces. J. Chem. Phys. 2005, 122, 104711. (20) Wang, J.; Kalinichev, A. G.; Kirkpatrick, R. J. Structure and Decompression Melting of a Novel, High-Pressure Nanoconfined 2-D Ice. J. Phys. Chem. B 2005, 109, 14308−14313. (21) Politano, A.; Marino, A. R.; Formoso, V.; Chiarello, G. Hydrogen Bonding at the Water/Quasi-Freestanding Graphene Interface. Carbon 2011, 49, 5180−5184. (22) Jinesh, K.; Frenken, J. Experimental Evidence for Ice Formation at Room Temperature. Phys. Rev. Lett. 2008, 101, 036101. (23) Kirov, M. V. New Two-Dimensional Ice Models. J. Stat. Phys. 2012, 149, 865−877. (24) Kirov, M. V. Residual Entropy of Ice Nanotubes and Ice Layers. Phys. A (Amsterdam, Neth.) 2013, 392, 680−688. (25) Anick, D. J. Atypical Water Lattices and their Possible Relevance to the Amorphous Ices: A Density Functional Study. AIP Adv. 2013, 3, 042119. (26) Baskin, Y.; Meyer, L. Lattice Constants of Graphite at Low Temperatures. Phys. Rev. 1955, 100, 544−544. (27) Haering, R. R. Band Structure of Rhombohedral Graphite. Can. J. Phys. 1958, 36, 352−362. (28) Leenaerts, O.; Partoens, B.; Peeters, F. Adsorption of H2O, NH3, CO, NO2, and NO on Graphene: A First-Principles Study. Phys. Rev. B 2008, 77, 125416. (29) Leenaerts, O.; Partoens, B.; Peeters, F. Water on Graphene: Hydrophobicity and Dipole Moment Using Density Functional Theory. Phys. Rev. B 2009, 79, 235440. (30) Ambrosetti, A.; Ancilotto, F.; Silvestrelli, P. L. Van der WaalsCorrected Ab Initio Study of Water Ice−Graphite Interaction. J. Phys. Chem. C 2013, 117, 321−325. (31) Dovesi, R.; Orlando, R.; Civalleri, B.; et al. CRYSTAL: A Computational Tool for the Ab Initio Study of the Electronic Properties of Crystals. Z. Kristallogr. 2005, 220, 571−573. (32) PQS Parallel Quantum Solutions, 2013 Green Acres Road, Fayetteville, AR 72703, USA. (33) Becke, A. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098−3100.

not a classic hydrophobe and also underscores that water can use a primarily vdW-based interaction to bond to a substrate. We saw that proton order matters at least somewhat, especially the T vs Z distinction. The very low barriers to shear motion predict easy gliding of HBL along graphene and may presage loss of templating for larger systems. For the ice-1h monolayer, there is a slight preference of “toward” over “away” for the nonH-bonded O−H’s, but the minimum energy occurs well short of 100%, with 55% toward being the ratio that we obtained. We saw how graphene’s cell dimensions do adapt a little to the templated H2O lattice, with the most usual situation being that more layers of water lattice induce greater graphene accommodation. The OHC motif occurs too consistently in the optimized geometries to be due to chance, but ultimately the difference it makes to graphene−HBL interactions may be minor. Our results support the preference for Gra+HBL over Gra +2L1h, as demonstrated by Kimmel et al., and suggest that once the ice is warm enough to rearrange H-bonds, the question of local ice thickness becomes a prominent issue for making accurate comparisons. Although single lattice calculations show a preference for Gra+4L1h over Gra+2HBL, the possibility of a two-lattice system and the various uncertainties about BLYP-D mean that the door is not entirely closed on Gra +2HBL. Lastly, we have lifted up Gra+nHBL as a structure worth looking for in confined geometries, via either experiment or modeling.

■

ASSOCIATED CONTENT

S Supporting Information *

Table comparing E0, E0 + cZPE, and cΔG298 for TAP computed with just the primitive cell and with various supercells. 2-D crystal coordinates for all the ice and graphene−ice structures listed in Tables 1 and 2. Figure 5 (illustrating four isomers of 2TAP) and Figure 6 (illustrating all-toward and all-away ice-1h on graphene). This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Notes

The author declares no competing financial interest.

■

REFERENCES

(1) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; et al. Electric Field Effect in Atomically Thin Carbon Films. Science 2004, 306, 666−669. (2) Novoselov, K. S.; Jiang, D.; Schedin, F.; et al. Two-Dimensional Atomic Crystals. Proc. Natl. Acad. Sci. U. S. A. 2005, 102, 10451− 10453. (3) Soldano, C.; Mahmood, A.; Dujardin, E. Production, Properties and Potential of Graphene. Carbon 2010, 48, 2127−2150. (4) Warner, J. H.; Schäffel, F.; Bachmatiuk, A.; Rummeli, M. H. Graphene: Fundamentals and Emergent Applications; Elsevier B.V.: Amsterdam, Netherlands, 2013. (5) Singh, V.; Joung, D.; Zhai, L.; et al. Graphene Based Materials: Past, Present and Future. Prog. Mater. Sci. 2011, 56, 1178−1271. (6) Liu, X.; Wang, M.; Zhang, S.; Pan, B. Application Potential of Carbon Nanotubes in Water Treatment: A Review. J. Environ. Sci. 2013, 25, 1263−1280. (7) Feller, D.; Jordan, K. D. Estimating the Strength of the Water/ Single-Layer Graphite Interaction. J. Phys. Chem. A 2000, 104, 9971− 9975. (8) Bermudez, V. M.; Robinson, J. T. Effects of Molecular Adsorption on the Electronic Structure of Single-Layer Graphene. Langmuir 2011, 27, 11026−11036. H

dx.doi.org/10.1021/jp500360n | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

(34) Grimme, S. Accurate Description of van der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (35) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (36) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab Initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (37) Tkatchenko, A.; Scheffler, M. Accurate Molecular van der Waals Interactions from Ground-State Electron Density and Free-Atom Reference Data. Phys. Rev. Lett. 2009, 102, 073005. (38) Von Lilienfeld, O. A.; Tkatchenko, A. Two- and Three-Body Interatomic Dispersion Energy Contributions to Binding in Molecules and Solids. J. Chem. Phys. 2010, 132, 234109. (39) Tkatchenko, A.; DiStasio, R. A.; Car, R.; Scheffler, M. Accurate and Efficient Method for Many-Body van der Waals Interactions. Phys. Rev. Lett. 2012, 108, 236402. (40) Reilly, A. M.; Tkatchenko, A. Seamless and Accurate Modeling of Organic Molecular Materials. J. Phys. Chem. Lett. 2013, 4, 1028− 1033. (41) Schäfer, A.; Huber, C.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets of Triple Zeta Valence Quality for Atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−5835. (42) Santra, B.; Michaelides, A.; Fuchs, M.; et al. On the Accuracy of Density-Functional Theory Exchange-Correlation Functionals for H Bonds in Small Water Clusters. II. The Water Hexamer and van der Waals Interactions. J. Chem. Phys. 2008, 129, 194111. (43) Tonigold, K.; Gross, A. Dispersive Interactions in Water Bilayers at Metallic Surfaces: A Comparison of the PBE and RPBE Functional Including Semiempirical Dispersion Corrections. J. Comput. Chem. 2012, 33, 695−701. (44) Yoo, S.; Xantheas, S. S. Communication: The Effect of Dispersion Corrections on the Melting Temperature of Liquid Water. J. Chem. Phys. 2011, 134, 121105. (45) Verdaguer, A.; Segura, J. J.; López-Mir, L.; Sauthier, G.; Fraxedas, J. Communication: Growing Room Temperature Ice with Graphene. J. Chem. Phys. 2013, 138, 121101. (46) Zheng, Y.; Su, C.; Lu, J.; Loh, K. P. Room-Temperature Ice Growth on Graphite Seeded by Nano-Graphene Oxide. Angew. Chem., Int. Ed. Engl. 2013, 52, 8708−8712. (47) http://pubs.acs.org. (48) Bayro-Corrochano, E.; Scheuermann, G. Geometric Algebra Computing: In Engineering and Computer Science; Springer-Verlag: London, U.K., 2010. (49) Rana, M. K.; Chandra, A. Ab Initio and Classical Molecular Dynamics Studies of the Structural and Dynamical Behavior of Water Near a Hydrophobic Graphene Sheet. J. Chem. Phys. 2013, 138, 204702. (50) Smith, D. E.; Dang, L. X. Computer Simulations of Cesium− Water Clusters: Do Ion−Water Clusters Form Gas-Phase Clathrates? J. Chem. Phys. 1994, 101, 7873−7881. (51) Anick, D. J. Topology-Energy Relationships and Lowest Energy Configurations for Pentagonal Dodecahedral (H2O)20X Clusters, X = Empty, H2O, NH3, H3O+: The Importance of O-Topology. J. Chem. Phys. 2010, 132, 164311. (52) Hirsch, T. K.; Ojamäe, L. Quantum-Chemical and Force-Field Investigations of Ice Ih: Computation of Proton-Ordered Structures and Prediction of Their Lattice Energies. J. Phys. Chem. B 2004, 108, 15856−15864. (53) Wei, S.; Shi, Z.; Castleman, A. W. Mixed Cluster Ions as a Structure Probe: Experimental Evidence for Clathrate Structure of (H2O)20 and (H2O)21H+. J. Chem. Phys. 1991, 94, 3268−3270. (54) Mishima, O. Reversible First-Order Transition Between Two H2O Amorphs at ∼0.2 GPa and ∼135 K. J. Chem. Phys. 1994, 100, 5910−5912. (55) Mundy, C. J. Personal communication, 2013.

I

dx.doi.org/10.1021/jp500360n | J. Phys. Chem. A XXXX, XXX, XXX−XXX