8412
J. Phys. Chem. C 2009, 113, 8412–8419
Structure and Stability of the Water-Graphite Complexes Miroslav Rubesˇ,† Petr Nachtigall,†,‡ Jirˇ´ı Vondra´sˇek,† and Ota Bludsky´*,† Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic and Center for Biomolecules and Complex Molecular Systems, FlemingoVo na´m. 2, CZ-16610 Prague, Czech Republic, and Department of Physical and Macromolecular Chemistry, Faculty of Science, Charles UniVersity in Prague, AlbertoV 6, 128 43 Praha 2, Czech Republic ReceiVed: February 16, 2009; ReVised Manuscript ReceiVed: April 6, 2009
The interaction of the water molecule with benzene, polycyclic aromatic hydrocarbons, graphene, and graphite is investigated at the density-functional/coupled-cluster (DFT/CC) level of theory. The accuracy of the DFT/ CC method is first demonstrated by a comparison of the various interaction energies on the potential energy surface of water-benzene, water-naphthalene, and water-anthracene complexes with the data calculated at the coupled-cluster level at the basis set limit. The potential energy surface of water-graphene and water-graphite is relatively flat with diffusion barriers of about 1 kJ/mol. The structure with both hydrogen atoms of water pointing toward the graphene plane (denoted as a circumflex structure) above the center of the six-member ring is the global minimum characterized with an electronic interaction energy of -13 and -15 kJ/mol for graphene and graphite, respectively. The OH · · · π complexes (with one OH pointing toward the surface and the other OH being oriented along the surface) are up to 2 kJ/mol less stable than the circumflex complexes of water on graphene/graphite, depending on the position of the oxygen atom. 1. Introduction The interaction of water molecules with various graphitic materials, including graphite, graphene, carbon nanotubes, fullerenes, and polycyclic aromatic hydrocarbons (PAHs), has recently attracted a great deal of interest among computational chemists.1-12 The motivation for this is relatively straightforward and can be characterized by the following attributes: (i) an understanding of the water interaction with carbon nanotubes or fullerenes is important for potential applications of these materials, (ii) the water-graphite and water-graphene represent model systems for the investigation of the interface between water and the hydrophobic surface, and (iii) the interaction of water with hydrophobic surfaces is relevant to the investigation of many biological systems. In biodisciplines, the interaction of water with hydrophobic surfaces is fundamental and interconnected with the hydrophobic effect demonstrated by structuring a water in the vicinity of nonpolar molecules.13 A large number of studies come from theory analyzing the structure of water at hydrophobic surfaces14,15 or during the hydrophobic collapse of a protein.16 Most of the theoretical works utilize a molecular dynamics simulation with approximate water models which cannot properly describe the subtle energy balance necessary for an accurate disclosure of the geometry arrangement. While the water interaction with graphite (or a related material) seems to be a simple issue, it represents a major challenge for both experimental and theoretical investigation. The fact that the water-water interaction is stronger than the water-graphite interaction is the main source of the complications on the experimental side. There are only a few relevant experimental studies addressing the interaction energy value of a single water molecule with the graphite surface. The * To whom correspondence should be addressed. E-mail: ota.bludsky@ uochb.cas.cz. Phone: +420-220410324. Fax: +420-220410320. † Academy of Sciences of the Czech Republic and Center for Biomolecules and Complex Molecular Systems. ‡ Charles University in Prague.
upper boundary for the water-graphite interaction comes from microcalorimetry, which yields a differential heat of adsorption at a zero-coverage limit of 19 kJ/mol;17 it is quite likely that the formation of water clusters on the surface results in excessively large interaction energy for an isolated water molecule. Other experimental values reported in the literature are the results of model calculations (employing empirical force fields) carried out to follow the experimentally observed property. Interaction energies in the range from -6.3 to -9.4 kJ/mol were obtained from the simulations1 aimed to reproduce the experimentally observed contact angles of water droplets on graphite.18,19 Obviously, there is a great demand for reliable model potentials describing the water interaction with graphite and similar materials. This in turn requires reliable data that could be used for the design and benchmarking of the model potentials. Highly accurate ab initio data, such as coupled clusters at the complete basis set limit (CCSD(T)/CBS), are currently available only for the water-benzene complex.20-22 The CCSD(T) results obtained with a medium-size basis set were reported for the water-naphthalene complex.3 However, the CCSD(T) calculations are computationally too expensive to be applied for larger PAH systems. Even if such calculations were available, the extrapolation of the interaction energies calculated for small PAHs to an infinite graphene (graphite) surface would represent a major challenge. The question of the PAH size that could represent graphene consistently has not been definitively answered. In this respect, a reliable and computationally tractable method must be used. The calculations at the MP2 level and at the density-functional theory (DFT) level have already been reported; while calculations on large PAHs were reported at the MP2 level the periodic calculations for graphene were reported at the DFT level. The accuracy of both of these methods for water-graphene interactions is yet to be established. The most extensive MP2 study of the water interaction with PAHs was carried out by Feller and Jordan.4 On the basis of
10.1021/jp901410m CCC: $40.75 2009 American Chemical Society Published on Web 04/20/2009
Structure and Stability of the Water-Graphite Complexes
J. Phys. Chem. C, Vol. 113, No. 19, 2009 8413
Figure 1. The definition of the reference set used for the DFT/CC calculations of the water-graphite interaction.
the extrapolation of the interaction energies of water complexes with benzene, coronene (C24H12), circumcoronene (C54H18), and dicircumcoronene (C96H24), they concluded that the watergraphene interaction amounts to -24 kJ/mol. This interaction energy appears to be too large in view of the experimental results as well as numerous simulations and has also been questioned by the authors themselves.5 Sudiarta and Geldart estimated the water-graphene interaction to be about -10 kJ/mol based on the extrapolation of the interaction energies of water-benzene, water-coronene, and water-circumcoronene at the MP2/631+G(0.25) level of theory; and similar interaction energy has been obtained from the extrapolation of the perfluorinated analogues.7 The DFT calculations also have been applied for the description of the water interaction with various PAHs. The traditional GGA-type exchange-correlation functionals cannot be used to describe weak intermolecular interactions (the water-graphite interaction is only about -4 kJ/mol).8,9 However, the extrapolation from PAHs to graphene/graphite is not necessary, because these methods can be applied to periodic systems. Provided that DFT can be accurately extended or corrected for weakly interacting complexes, it can be advantageously used for the description of periodic surfaces, such as graphene, graphite, or carbon nanotubes, employing the planewave basis set. There has been great progress recently in extending the DFT formalism to account for the weak intermolecular interactions. The interaction energy of the waterbenzene complex of -12 kJ/mol was found with the vdW-DF functional,23 slightly below the CCSD(T)/CBS results. A density-functional tight-binding method augmented with an empirical van der Waals correction (DFTB-D) was also used for the description of the water interaction with benzene and various PAHs (up to dicircumcoronene).10 The DFT-D approach has been evaluated for the water-benzene complex considering different functionals, damping functions, and C6 coefficients;24 it was concluded that the B3LYP functional with dispersion correction is the most accurate one (the deviation from the reference level of theory not exceeding 2 kJ/mol). Note, however, that the reference level of theory used in ref 24, MP2/ 6-311G(2d,2p), yields a result deviating from the accurate CCSD(T)/CBS level (e.g., the stability of the parallel complex that is unbound at the MP2/6-311G(2d,2p) level). A recently developed density-functional theory/coupledcluster method (DFT/CC)25 has been applied here for the description of the water interaction with PAHs, graphene, and graphite. The basic concept of the DFT/CC method is the following: the DFT interaction energies of weakly bound systems can be corrected for the CCSD(T)/CBS accuracy on a suitable reference set of molecules. The correction is evaluated for one-dimensional curves on the potential energy surface (PES), without any assumption about the functional form of this correction. The correction is then used for larger systems of a similar type (including periodic ones) assuming a good transferability of the correction. This method has been recently used in the investigation of various weakly interacting hydro-
carbons in the gas and solid phases, for which a very good accuracy of the method was reported.26-28 The extension of the method to systems containing polar molecules is reported here. The DFT/CC method is described with the reference set (water-benzene and water-molecular hydrogen) being defined first (Section 2). The accuracy of the DFT/CC method is analyzed in Section 3.1 and then used in the investigation of water interacting with various larger PAHs and with graphene and graphite (Sections 3.2 and 3.3, respectively). 2. Methods The DFT/CC correction scheme has been developed for an accurate description of weakly bound molecular systems.25,27 The underlying assumption of the method is the pairwise representability of the DFT error, ∆EDFT/CC, in terms of the interatomic distances, Rij,
∆EDFT⁄CCSD(T) )
∑ εij(Rij) ij
where εij are as yet unknown correction functions. In the DFT/ CC methodology, ∆EDFT/CC is defined as the difference between the CCSD(T)/CBS and DFT interaction energies. The pairwise representability assumption holds remarkably well for weakly bound complexes, as has been recently shown for hydrocarbon systems.25-28 The correction functions are obtained from the potential energy functions calculated at the CCSD(T) and DFT levels (employing the Perdew-Burke-Ernzerhof, PBE,29 exchange-correlation functional) for a suitable reference set of molecules. In this work, the DFT/CC correction functions were obtained from the reference set consisting of H2 · · · H2O and H2O · · · C6H6 (at the geometries depicted in Figure 1) by means of the Reciprocal Power Reproducing Kernel Hilbert Space Interpolation (RP-RKHS).30 No a priori functional form of the correction functions is assumed (cf. DFT-D) except for the asymptotic behavior given by the RKHS kernel (R-6 + R-8). A more detailed description of the DFT/CC method can be found in refs 25 and 27. The correction functions εij are reported in the Supporting Information; additional details about the εij evaluation are also reported in the Supporting Information. The accuracy of the DFT/CC method relies on the precise (CCSD(T)/CBS) calculations for the reference systems (waterhydrogen and water-benzene). The performance of the DFT/ CC method is checked based on the comparison of the DFT/ CC results with the CCSD(T)/CBS results for larger complexes (water-naphthalene and water-anthracene). Due to the computational demands the way the extrapolation on the CBS limit is done depends on the size of these systems. The Dunning’s correlation-consistent valence-X-ζ basis sets with polarization functions with X ) double, triple, quadruple or quintuple basis sets are denoted as VDZ, VTZ, VQZ, and V5Z, respectively.31 The corresponding augmented Dunning’s correlation-consistent basis sets are labeled as AVDZ, AVTZ, AVQZ, and AV5Z, respectively. The H2 · · · H2O complex was calculated at the CCSD(T)/AVTZ and CCSD(T)/AVQZ level,
8414
J. Phys. Chem. C, Vol. 113, No. 19, 2009
Rubesˇ et al.
and the CCSD(T)/CBS estimate was obtained from the simple correlation energy dependence on the cardinal number X (EX ) ECBS + AX -3), where the HF/AV5Z energy was taken as the CBS limit for the uncorrelated part. Water-benzene was calculated at the CCSD(T)/AVTZ level and the CCSD(T)/CBS was obtained from the following formula
CCSD(T)/CBS ) (CCSD(T)/AVTZ) + (MP2/CBS) (MP2/AVTZ) where the MP2/CBS estimate was obtained from the MP2/ AVQZ and MP2/AVTZ by means of the simple extrapolation formula mentioned above and HF/AV5Z* (AV5Z on water and carbon atoms and of V5Z on H atoms of a PAH) was used as the CBS limit for the uncorrelated part. The water-naphthalene and water-anthracene complexes were calculated at the CCSD(T)/AVDZ level and subsequent extrapolation to the CCDS(T)/ CBS limit was performed according to the formula
CCSD(T)/CBS ) (CCSD(T)/AVDZ) + (MP2/CBS) (MP2/AVDZ) where the MP2/CBS estimate was obtained in the same manner as for the water-benzene complex by using the HF/AV5Z** (V5Z on all H atoms) level for the uncorrelated part. The DFT calculations for the reference set were performed with the PBE functional and AVQZ basis set. The PBE calculations for larger water-PAH complexes were performed with the VTZ and AVTZ basis sets on hydrocarbon and water, respectively, denoted as AVTZ/VTZ. The resolution of identity,32 employing the auxiliary basis sets from refs 33 and 34, were used in all the calculations at the DFT and MP2 levels. To avoid the linear-dependency problems observed for MP2 calculations of coronene,4 the mixed basis consisting of AVTZ on six central C atoms and on water and of VTZ on other PAH atoms (labeled as (A)VTZ) was used. All the calculated interaction energies of the water-PAH complexes were corrected for the basis set superposition error by means of the counterpoise correction of Boys and Bernardi.35 Most of the calculations were performed within the frozen monomer approximation; the CCSD(T)/AVQZ geometries of H2 (rHH ) 0.7420 Å) and H2O (C2V) (rOH ) 0.9590 Å, R ) 104.4°) and the geometry of Gauss and Stanton36 for C6H6 (D6h) (rCC ) 1.3915 Å, rCH ) 1.0800 Å) were used. The naphthalene and anthracene monomers (D2h) were obtained from a geometry optimization at the PBE/AVQZ level. The geometries of the larger PAHs (coronene, circumcoronene, and dicircumcoronene) were obtained at the PBE/VTZ level. The plane-wave PBE calculations were carried out for the water-PAH (the circumflex structure) and water-graphene/ graphite complexes. The Kohn-Sham equations were solved through the projected augmented wave method of Blo¨chl.37 The plane-wave basis set with an energy cutoff of 800 eV and the PAW pseudopotential with ENMAX ) 700 eV were used. The plane-wave water-PAH calculations were performed by using a 20 × 20 × 20 Å3 simple cubic cell for benzene, naphthalene, and coronene and a 25 × 25 × 20 Å3 cell for circumcoronene. The calculations on graphene were performed with an interlayer separation of 20 Å and with the PBEoptimized C-C bond length (1.422 Å). To prevent water from interacting with its images in the periodic calculations, the supercell consisting of 72 carbon atoms was adopted with a lateral water-water distance of about 15 Å. The experimental value of the c parameter (6.67 Å) was used for the calculations
Figure 2. The water-PAH and water-graphene/graphite structures defined by the water orientations (I-IV) and the positions of the oxygen atom above the six-member carbon ring (a-c).
on graphite. The Γ-point calculations were found to be sufficient forobtainingfullyconvergedinteractionenergiesofwater-graphene/ graphite. The following program suites were used: Turbomole 5.9,38 Gaussian03,39 and Molpro 2006.1.40 The DFT/CC optimizations were performed through the link External in Gaussian03. The periodic DFT calculations were performed with the Vienna ab initio Simulation Package (VASP).41 3. Results The water-PAH and water-graphene/graphite structures investigated in this study are depicted in Figure 2. The structures are described by water orientations (I, circumflex; II, OH · · · π; III, coplanar; and IV, O-down structures) and the positions of the oxygen atom over the six-member carbon ring (a, over the center of the ring; b, over the middle of the C-C bond; and c, over the atom). The water C2 axis was fixed parallel to the PAH normal in structures I and IV; the O-H bond was fixed along the PAH normal in structure II. Unless stated otherwise, the geometry optimizations were performed within the constraints defined above. Note that some of the structures defined above are not necessarily stationary points on the PES (e.g., the OH · · · π structure IIa is not a stationary point on the waterbenzene PES); however, it is useful to compare the stabilities of such constrained structures as they depend on the size of the PAH. Considering the fact that the frozen monomer geometries were used, it follows that the structures were obtained by a onedimensional optimization. 3.1. The Reliability of the DFT/CC Method. The equilibrium intermolecular distances and interaction energies of the water complexes with benzene, naphthalene, and coronene obtained at the DFT/CC, MP2/AVTZ, and CCSD(T)/CBS levels are listed in Table 1. The transferability of the pairwise correction functions is a crucial issue for the reliability of all empirically corrected DFT methods. The choice of reference system is thus the main factor influencing the quality of the DFT/CC predictions. In this work, the correction functions were constructed from the CCSD(T) and DFT potential energy curves calculated for selected structures of the water-H2 and water-benzene complexes (Ia and IVa, see Figure 1). It is then obvious that the DFT/CC method yields interaction energies for the water-benzene complex that are in very good agreement with the CCSD(T)/CBS level of theory (see Table 1). Note, however, that good agreement between the CCSD(T) and DFT/ CC was found also for the water-benzene complexes of types II and III, which were not included in the reference set (Table 1). To test the transferability of the correction functions, coupledcluster calculations onto the water-naphthalene complexes were
Structure and Stability of the Water-Graphite Complexes TABLE 1: The Equilibrium Intermolecular Distances (in Å) and Interaction Energies (in kJ/mol) of Benzene, Naphthalene, and Coronene Complexes with Water Calculated at the MP2/AVTZ, CCSD(T)/CBS, and DFT/ CC(AVQZ) Levels MP2/AVTZb structure
a
Re
Eint
CCSD(T)/CBS Re
Eint
DFT/CC(AVQZ) Re
Eint
Ia IIa IIIa IVa
benzene · · · H2O 3.32 -13.2 3.32 -13.2 3.37 -13.2 3.38 -12.9 3.43 -2.2 3.43 -2.0 unbound unbound
Ia Ib IIa IIb IIIb IVb
3.28 3.30 3.38 3.39 3.23 3.21
naphthalene · · · H2O -13.9 3.29 -13.3 -15.3 3.31 -14.7 -12.3 3.40 -11.6 -13.2 3.42 -12.4 -6.0 3.27 -5.2 -2.2 3.27 -1.0
3.28 3.33 3.41 3.41 3.25 3.27
-13.6 -14.3 -11.3 -12.2 -5.2 -0.7
Ia IIa IIb IIc IIIa IVa
3.26 3.32 3.36 3.36 3.11 2.98
coronene · · · H2O -16.1 -13.7 -14.3 -14.5 -10.3 -7.2
3.27 3.37 3.38 3.37 3.12 3.04
-14.8 -11.8 -12.8 -13.0 -8.8 -4.8
3.32 -13.2 3.40 -12.3 3.41 -2.2 unbound
a
The structure types are defined in Figure 2. b The (A)VTZ basis was used for MP2 calculations on water-coronene.
performed and the results compared to those obtained at the DFT/CC level. The potential energy curves (Ib, IIa, IIb, IIIb, and IVb structures) calculated at the DFT/AVQZ, CCSD(T)/ CBS, and DFT/CC levels of theory are depicted in Figure 3. At the equilibrium intermolecular separation, Re, the maximum energy difference between the CCSD(T) and DFT/CC values is about 0.4 kJ/mol and the values of Re obtained at both levels of theory are within 0.02 Å. Note that the DFT/CC method predicts correct interaction energies along the entire potential energy curves. It can be concluded that the transferability of the DFT/CC corrections is good and that the PES for both the water-benzene and water-naphthalene complexes can be reliably reproduced by the DFT/CC method (to within a few tenths of a kilojoule per mole). This conclusion is further supported by the CCSD(T)/CBS and DFT/CC calculations for the Ia structure of the water-anthracene complex. The CCSD(T)/ CBS value of -13.5 kJ/mol differs by 0.4 kJ/mol from the DFT/ CC estimate of -13.9 kJ/mol. While the global PES of the relatively small water-PAH complexes (benzene or naphthalene in particular) can be currently explored at the CSSD(T)/CBS level, the use of the CCSD(T) method for larger PAHs, like coronene, is computationally too expensive. Therefore, it is meaningful to compare the DFT/CC (and CCSD(T) where available) results with another correlated ab initio method, such as MP2. It has been reported previously that the MP2 method performs reasonably well for the water-benzene complex.20 An excellent agreement between the MP2/AVTZ and CCSD(T)/CBS results was found for the water-benzene complexes reported in Table 1. With respect to the CCSD(T)/CBS level, the MP2/AVTZ gives interaction energies that are slightly overestimated. A similar trend is also apparent for the water-naphthalene complex, but the differences between the CCSD(T)/CBS and MP2/AVTZ interaction energies become larger with the greatest deviation found for the IVb complex, the interaction energy of which is overestimated by 1.2 kJ/mol at the MP2/AVTZ level (contrary
J. Phys. Chem. C, Vol. 113, No. 19, 2009 8415 to the DFT/CC interaction energies, which are within 0.4 kJ/ mol of the CCSD(T)/CBS values). The MP2/AVTZ interaction energies for water-naphthalene are always larger than the DFT/ CC results; the difference is in the range from 0.3 to 1.5 kJ/ mol. This difference becomes larger (1.3-2.4 kJ/mol) for the water-coronene complexes. Assuming that the CCSD(T)/CBS results represent the most reliable estimates of water-PAH interaction currently available, it can be concluded that the DFT/ CC method provides reliable results for these complexes while the accuracy of the MP2 method deteriorates with increasing PAH size. 3.2. The Water-PAH Complexes. The effect of the PAH size on the interaction energy with the water molecule was investigated at the DFT/CC level for the circumflex-structure type (Table 2). The motivation for this selection is the following: (i) the circumflex structure is more stable than the OH · · · π structure for the water-graphene complex (Table 3), (ii) the OH · · · π complex (that is a global minimum on the water-benzene PES) is not a minimum on the water-graphene PES, and (iii) the increased stability of the circumflex structure with respect to the OH · · · π structure is already apparent for the water-coronene system (Table 1). The water-PAH interaction energies and equilibrium distances calculated at the DFT/CC level for benzene, naphthalene, anthracene, coronene, circumcoronene, dicircumcoronene (37 fused benzene rings), and infinite (periodic) graphene are reported in Table 2; calculations employing the AVQZ (only up to the size of coronene, because the linear dependence problem means that convergence was not achieved for larger systems), AVTZ/VTZ, and PW basis sets are reported. The results obtained with the AVQZ and PW basis sets are in excellent agreement. Also the results obtained with the smaller AVTZ/VTZ basis set are in very good agreement with those obtained with the PW basis set (the largest difference being 0.3 kJ/mol). Very good agreement between the results obtained with the plane-wave and atom-centered basis sets allows the conclusion to be drawn that the results obtained at the DFT/ CC level employing periodic models have the same accuracy as those obtained for the finite-size PAHs. The interaction energy obtained at the PBE level (reported in parentheses in Table 2) monotonically decreases when going from water-benzene to water-graphene; correspondingly, the equilibrium distance increases from water-benzene to watergraphene. On the contrary, the interaction energy calculated at the DFT/CC level increases from water-benzene to watercoronene and then decreases with the increasing size of the PAH; note that the equilibrium distance monotonically decreases from water-benzene to water-graphene. This observation can be rationalized as follows: the intersystem correlation term that is missing at the DFT level (often referred to as dispersion) and that is represented by the εij correction function decays quickly with increasing water-graphene distance. The correction term for C atoms on the outer ring of circumcoronene (about 1 kJ/mol) is smaller than the change in the DFT interaction energy between coronene and circumcoronene (1.2 kJ/mol). Thus, the DFT/CC correction compensating mainly for the attractive dispersion component of the water-PAH interaction reverses the trend for the small PAH systems (benzene, naphthalene, and coronene) while for larger PAHs the interaction energy starts to decrease again as the DFT trend becomes dominant (Figure 4). 3.3. The Water-Graphene/Graphite System. The equilibrium distances and interaction energies for various complexes of water on graphene and graphite are given in Table 3. The DFT part of the water-graphite DFT/CC interaction energies is taken from the calculations on water-graphene (i.e., at the
8416
J. Phys. Chem. C, Vol. 113, No. 19, 2009
Rubesˇ et al.
Figure 3. A comparison of the DFT, CCSD(T), and DFT/CC potential energy curves for water-naphthalene complexes as a function of intermolecular separations (with R being the distance between the oxygen atom and the naphthalene plane).
TABLE 2: The DFT/CC Equilibrium Distances (in Å) and Interaction Energies (in kJ/mol) of the Circumflex Structure Calculated for Water-PAH Complexes of Increasing Sizea AVTZ/VTZc
AVQZ n
b
1 2 3 7 19 37 graphene
PW
Re
Eint
Re
Eint
Re
Eint
3.32 (3.50) 3.28 (3.50) 3.26 (3.51) 3.27 (3.62)
-13.2 (-7.7) -13.6 (-7.1) -13.9 (-6.3) -14.8 (-5.7)
3.32 (3.51) 3.29 (3.51) 3.26 (3.52) 3.28 (3.64) 3.28 (3.66) 3.27 (3.66)
-13.4 (-7.9) -13.8 (-7.3) -14.0 (-6.5) -15.0 (-5.9) -14.7 (-4.6) -14.5 (-4.2)
3.32 (3.50) 3.30 (3.50) 3.27 (3.51) 3.28 (3.62) 3.27 (3.63)
-13.2 (-7.7) -13.7 (-7.2) -13.9 (-6.3) -14.7 (-5.6) -14.5 (-4.4)
3.26 (3.65)
-13.3 (-2.7)
a The DFT (PBE) results are given in parentheses. b The number of the PAH’s benzene rings. c The VTZ basis on a PAH and the AVTZ basis on water.
TABLE 3: The DFT/CC Equilibrium Intermolecular Distances (in Å) and Interaction Energies (in kJ/mol) of Water-Graphene and Water-Graphite Complexes graphene
graphite
structurea
Re
Eintb
Re
Eint
Ia Ib Ic IIa IIb IIc IIIa IVa
3.26 3.34 3.32 3.35 3.35 3.35 3.12 3.01
-13.3 (-1.1) -12.0 (-0.6) -12.4 (-0.9) -11.4 (-1.4) -11.9 (-2.0) -12.1 (-2.1) -10.0 (0.4) -8.7 (-0.1)
3.24 3.32 3.31 3.34 3.34 3.34 3.11 3.00
-15.0 -13.7 -14.1 -12.7 -13.3 -13.5 -10.8 -8.8
a The structure types are defined in Figure 2. parentheses are the PBE interaction energies.
b
The values in
DFT level, the interaction energies for graphene and graphite are taken to be the same). The error of this approximation, however, is very small (within 0.1 kJ/mol), because the DFT calculations on graphene and on graphite, represented by two graphene sheets, do not differ by more than 0.1 kJ/mol. The changes between the graphene and graphite PESs are rather subtle. The energy differences between the reported structures are slightly larger for graphite, but the shape of the
PES is preserved. The lowest energy arrangements of water are stabilized by forming OH · · · π bonds (I and II), whereas the O-down structures (IV) are the least favorable (-8.8 kJ/mol). The coplanar structure (III) with an interaction energy of -10.8 kJ/mol can be viewed as a transition between the circumflex (I) and O-down (IV) stationary structures. Among the circumflex (I) structures investigated on the water-graphene PES, the structure above the center of the ring (Ia) is the global minimum while the structure above the center of the C-C bond (Ib) is the transition-state structure. The PES around the circumflex structure above the atom (Ic) is rather flat; the Ic structure may not be the stationary point, but its energy is close to a transitionstate structure of the first order connecting two adjacent Ia minima. Apparently, the water-graphite PES is rather flat with low barriers (∼1 kJ/mol) to the lateral motion of water above the graphite surface. The most stable adsorption complex of water on graphite and graphene is the circumflex complex above the center of the sixmember ring (-15.0 and -13.3 kJ/mol, respectively). The OH · · · π complex over the center of the ring (IIa) found to be the global minimum on the water-benzene PES is not a minimum on the water-graphene/graphite PESs. The 1-D constrained optimization of the IIa complex yields an interaction
Structure and Stability of the Water-Graphite Complexes
J. Phys. Chem. C, Vol. 113, No. 19, 2009 8417
Figure 4. The interaction energies of the water-PAH and water-graphene complexes (circumflex structure) depending on the PAH size. The PBE, DFT/CC, and MP2 results are depicted in green, blue, and red, respectively; the CCSD(T)/CBS results are given as stars.
energy of -11.4 kJ/mol for water-graphene, less than that found for the corresponding complex between water and benzene. On the contrary, the circumflex structure on the water-benzene PES is less stable than any of the circumflex structures on the water-graphene PES. For the circumflex structure, the minimum is over the center of the carbon ring (Ia). Interestingly, the most stable site for the asymmetric OH · · · π structure is over the carbon atom (IIc). 4. Discussion Previous estimates of water-graphite interactions have been scattered in the range from -6 to -24 kJ/mol. The upper limit, -24 kJ/mol reported by Feller and Jordan,4 is based on an extrapolation of the MP2 interaction energies calculated for PAHs with increasing size. The lower limit includes the interaction energies determined indirectly from the experimental contact angles of water droplets on the graphite surface. A rather broad range of values, ranging from 42° to 86°, was reported;18,19 using MD simulation, it was shown that the contact angle changes significantly as a function of the water-carbon interaction energy.1 To recover the macroscopic contact angles of 86° and 50°, a water monomer binding energy of -6.3 and -9.4 kJ/mol, respectively, was determined as necessary. The numerous empirical potentials used in the MD studies of water-graphite interaction yield values of about -10 kJ/mol.5 To narrow the range of the estimated interaction energies for water-graphite, it is essential to assess the accuracy of the methods used for the system’s characterization. The complexes of water with small PAHs (benzene, naphthalene, anthracene) are the most suitable benchmark systems currently available as they can provide an assessment of the various theoretical approaches with respect to the most accurate and reliable CCSD(T)/CBS estimates of water-PAH binding. The accuracies of the methods for the equilibrium structure determination are discussed first. The geometries fully optimized at the CCSD(T)/CBS level are not available. One-dimensional curves on the PES were calculated (within the frozen monomer
approximation) for selected structures of water-benzene and water-naphthalene at the CCSD(T)/CBS level and are in very good agreement with the corresponding results obtained at the DFT/CC level (differences within 0.02 Å, Table 1 and Figure 3). The interaction energy of the global minimum structure of the water-benzene complex fully optimized at the DFT/CC level is -13.5 kJ/mol while the CCSD(T)/CBS at this geometry yields -13.9 kJ/mol. Note that at the global minimum optimized at the DFT/CC level the water-benzene complex has the OH · · · π structure with one O-H oriented toward the benzene ring. This is in good agreement with the previous MP2 and CCSD(T) optimizations20 and with the experimentally determined structure.42From the microwave spectra, it was shown that θ, defined as an angle between the benzene normal and the water C2 axis, is 20((15)°, indicating that even the barriers separating the minima on the PES for benzene (a circumflex structure) are rather low (the DFT/CC gave θ ) 25° and 22°, with and without the frozen monomer approximation, respectively).42,43 The binding energy of the water-benzene complex calculated with the DFT/ CC (-13.5 kJ/mol) is also in good agreement with the CCSD(T)/ CBS estimates published previously (-14.1 kJ/mol44 and -13.7 kJ/mol22) and with the available experimental results -14.4 kJ/ mol45 and -13.6 kJ/mol46 (considering the calculated ZPE20 correction). The performance of a computational method when passing from the water-benzene complex to more extended systems is of crucial importance for the accuracy of water-graphite calculations. The water complexes with benzene, naphthalene, and anthracene, for which the CCSD(T)/CBS interaction energies were calculated, represent a unique set for the assessment of methods and extrapolation accuracy (only the Ia structure, circumflex above the center of the ring, has been calculated at the CCSD(T)/CBS level for water-anthracene). In addition to good agreement between the DFT/CC and CCSD(T)/CBS interaction energies calculated for various water-benzene and water-naphthalene complexes (Table 1), it is important that the DFT/CC error does not depend significantly on the size of
8418
J. Phys. Chem. C, Vol. 113, No. 19, 2009
the PAH. For the Ia structure, the differences between the CCSD(T) and DFT/CC are 0.3 and 0.4 kJ/mol for naphthalene and anthracene complexes, respectively (Figure 4). The consistent behavior of the DFT/CC and CCSD(T)/CBS methods provides confidence in the water interaction energy with graphene and graphite (13.3 and 15.0 kJ/mol, respectively). Note that these values do not result from the extrapolation of the interaction energies calculated for the PAHs. In addition, the DFT/CC method within the periodic framework has been clearly established (Table 2). Thus, the accuracy of the adsorption enthalpies calculated for the water-graphite system depends only on the transferability of the εij correction functions; on the basis of the results obtained for water complexes with benzene, naphthalene, and anthracene, it can be concluded that the transferability of the εij correction functions is reasonably good. It is also useful to discuss the accuracy of the MP2 method, which is often used in the investigation of weak intermolecular interactions. The overestimating tendency of the MP2/AVTZ method is apparent already on the water-naphthalene and water-anthracene complexes. For circumflex (I) and OH · · · π (II) structures of water-naphthalene, the MP2 overestimates rather systematically by 0.6-0.7 kJ/mol. Moreover, the difference between the MP2 and CCSD(T) increases with the size of the PAH (see Figure 4). For the Ia structure of anthracene, the MP2 interaction energy of -14.6 kJ/mol is 1.1 kJ/mol larger than the corresponding CCSD(T)/CBS value (-13.5 kJ/mol). The same trend is observed for coronene, where the MP2 value of -16.1 kJ/mol is likely to be overestimated by 1.5-1.7 kJ/ mol for the circumflex structure of the water-coronene complex. The water-PAH and water-graphite complexes have been investigated at the MP2/CBS level of theory.4 The impressive work on extended water-PAH complexes as large as C96H24 predicts the extrapolated MP2 interaction energy for watergraphene of -24 ( 2 kJ/mol. This value differs significantly (by 11 kJ/mol) from the DFT/CC prediction of -13 kJ/mol for this system. There are two main sources of this discrepancy: (i) the excellent performance of the MP2 for small model systems like water-benzene complex deteriorates for larger PAHs (see above) and (ii) the extrapolation procedure of ref 4 is based on monotonic behavior toward the infinite water-PAH (graphene) system. While the first effect leads to only a small overestimation of the MP2 binding within a few kilojoules per mole (1.3 kJ/mol for coronene, Figure 4) and hence cannot account for the roughly 11 kJ/mol difference between the MP2 and DFT/CC estimates for water-graphene, the second effect leads to a much more serious discrepancy. As discussed in Section 3.2, the DFT/CC predicts a decrease in the interaction energy for very large PAH systems in contrast to the reported MP2 tendency to increase the water-PAH binding monotonically. Unfortunately, there is no computationally tractable and reliable method to resolve this contradiction conclusively. As can be seen from Figure 4, the smallest system to check the asymptotic behavior of the MP2 and DFT/CC with respect to the size of PAH systems is circumcoronene (C54H18), which is already too large for highly correlated ab initio calculations such as the CCSD(T). Thus, the only arguments to support the DFT/ CC results are indirect. Note, however, that our estimate of the water-graphene interaction (-13.3 kJ/mol) is in good agreement with the results of the independent study of Jenness and Jordan.47 On the basis of the analysis of the DFT-SAPT investigation of water interaction with PAHs, they concluded that the water-graphene interaction for an OH · · · π complex above the center of the ring was -10.4 kJ/mol. This binding energy is within 1 kJ/mol of our results reported for the
Rubesˇ et al. corresponding (IIa) structure (-11.4 kJ/mol, Table 3). Part of the discrepancy is likely due to the differences in geometries used in both studies. The convergence of the interaction energies with the PAH size appears to be rather slow at the PBE and DFT/CC levels (Figure 4). The interaction energy of water with graphene is close to that of water with benzene. However, the PESs of water-benzene and water-graphene are qualitatively different. In light of our results reported in Tables 1 and 2, the modeling of water-graphite interaction by using medium-size PAHs seems to be questionable. The full periodic DFT calculations employing a plane-wave basis set are relatively straightforward, avoiding thus the uncertainty in the extrapolations. The use of the currently implemented local functionals, however, results in the recovery of only a small fraction of the total interaction energy, which must therefore be corrected for the missing part of the correlation energy. The interaction energies calculated for water-graphene at the DFT and DFT/CC levels (Table 3) show that the DFT not only underestimates the interaction by more than 10 kJ/ mol but also gives qualitatively incorrect results for the relative stabilities of the various complex types (e.g., the most stable circumflex complex Ia is about 1.2 kJ/mol more stable than the most stable OH · · · π complex (IIc) at the DFT/CC level while the IIc complex is 1 kJ/mol more stable than the Ia complex at the DFT level). It follows that the results obtained at the DFT level employing only GGA-type functionals11,12 should be taken with caution. A failure of the DFT calculation is also apparent from a comparison of the results for the singlelayer (graphene) and two-layer (graphite) models. For the circumflex structure (Ia), the addition of the second layer results in a tiny increase of the interaction energy by less than 0.1 kJ/ mol while the DFT/CC method results in an increase of 1.6-1.8 kJ/mol when going from the single-layer to multilayer model. When properly tested, the DFT/CC method has been shown to provide a reliable framework for obtaining very accurate estimates of the equilibrium geometries and interaction energies of weakly bound molecular systems.25-28 In the case of water interacting with the graphitic surface, it is estimated that the DFT/CC overestimates the interaction, but this error should not exceed 0.5 kJ/mol. The calculated interaction energy of water is about -13 and -15 kJ/mol for graphene and graphite, respectively. It is concluded that D0 for the water interaction with graphene and graphite is -9 and -11 kJ/mol, respectively (applying the ZPE estimate from ref 20). 5. Conclusions The interaction of water with polycyclic aromatic hydrocarbons, graphene, and graphite was investigated at the DFT/CC level of theory. The accuracy of the DFT/CC method was examined based on a comparison with the CCSD(T)/CBS results. Very good agreement was found for the interaction energies of water with benzene, naphthalene, and anthracene; in addition, the difference between the CCSD(T)/CBS and DFT/ CC results does not depend on the PAH size. This agreement allows for the use of the DFT/CC method in the investigation of the water interaction with larger PAHs and graphite; the estimated error for the calculated interaction energies was 0.5 kJ/mol. The calculations on water-graphite show that the circumflex structure above the center of the six-membered ring, characterized by an interaction energy of -15.1 kJ/mol, is a global minimum structure on the water-graphite PES. Both O-H bonds are equivalent in the circumflex structure and are oriented
Structure and Stability of the Water-Graphite Complexes toward the graphite surface. The O-down structure is energetically the least favorable (Eint ) -8.8 kJ/mol). The coplanar structure with an interaction energy of -10.8 kJ/mol can be viewed as a transition between the circumflex and O-down stationary structures on the global water-graphite PES. The OH · · · π structure with one OH pointing toward the surface is less stable (Eint ) -13.5 kJ/mol) than the circumflex structure, whereas the OH · · · π structure above the C atom is more stable than the one above the center of the six-member ring. These results are qualitatively different from those found for the water-benzene complex, where the OH · · · π complex above the center of the six-member ring was found to be the global minimum. The results reported here, along with those reported by Jenness and Jordan,47 constitute a benchmark for other methods that should be used in future studies. The results for the water-graphene complexes are very similar to those for water-graphite; the subsurface graphene sheets provide additional stabilization of the water-graphene complexes (0.1-1.7 kJ/mol stabilization depending on the structure type). Both PESs are rather flat with low barriers (∼1 kJ/mol) to the lateral motion of water over the surface. Acknowledgment. This research has been supported by grant No. IAA400550613 of the Grant Agency of the ASCR. Work at the Institute of Organic Chemistry and Biochemistry has been supported by the Ministry of Education, Youth and Sports (MEYS) grant No. LC512 and research project No. Z4 055 0506. Work at Charles University has been supported by the MEYS research project No. 0021620857. We would like to thank Ken Jordan for sending us a preprint of his paper on water-graphene interaction. Supporting Information Available: The DFT/CC correction functions along with additional details about their evaluation. This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Werder, T.; Walther, J. H.; Jaffe, R. L.; Halicioglu, T.; Koumoutsakos, P. J. Phys. Chem. B 2003, 107, 1345. (2) Werder, T.; Walther, J. H.; Jaffe, R. L.; Koumoutsakos, P. Nanotech 2003, 3, 2003–546. (3) Cabaleiro-Lago, E. A.; Rodriguez-Otero, J.; Pena-Gallego, A. J. Phys. Chem. A 2008, 112, 6344. (4) Feller, D.; Jordan, K. D. J. Phys. Chem. A 2000, 104, 9971. (5) Karapetian, K.; Jordan, K. D. Water in Confining Geometries; Buch, V., Devlin, J. P., Eds.; Springer-Verlag: Berlin, Germany, 2003; p 139. (6) Pertsin, A.; Grunze, M. J. Phys. Chem. B 2004, 108, 1357. (7) Sudiarta, I. W.; Geldart, D. J. W. J. Phys. Chem. A 2006, 110, 10501. (8) Leenaerts, O.; Partoens, B.; Peeters, F. M. Phys. ReV. B 2008, 77. (9) Wehling, T. O.; Lichtenstein, A. I.; Katsnelson, M. I. Appl. Phys. Lett. 2008, 93. (10) Lin, C. S.; Zhang, R. Q.; Lee, S. T.; Elstner, M.; Frauenheim, T.; Wan, L. J. J. Phys. Chem. B 2005, 109, 14183. (11) Sanfelix, P. C.; Holloway, S.; Kolasinski, K. W.; Darling, G. R. Surf. Sci. 2003, 532, 166. (12) Cicero, G.; Grossman, J. C.; Schwegler, E.; Gygi, F.; Galli, G. J. Am. Chem. Soc. 2008, 130, 1871. (13) Scatena, L. F.; Brown, M. G.; Richmond, G. L. Science 2001, 292, 908. (14) Raschke, T. M.; Levitt, M. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 6777.
J. Phys. Chem. C, Vol. 113, No. 19, 2009 8419 (15) Raschke, T. M.; Tsai, J.; Levitt, M. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 5965. (16) Zhou, R. H.; Huang, X. H.; Margulis, C. J.; Berne, B. J. Science 2004, 305, 1605. (17) Avgul, N. N.; Kieslev, A. V. Chemistry and Physics of Carbon; Dekker: New York, 1970. (18) Fowkes, F. M.; Harkins, W. D. J. Am. Chem. Soc. 1940, 62, 3377. (19) Schrader, M. E. J. Phys. Chem. 1980, 84, 2774. (20) Feller, D. J. Phys. Chem. A 1999, 103, 7558. (21) Tsuzuki, S.; Honda, K.; Uchimaru, T.; Mikami, M.; Tanabe, K. J. Am. Chem. Soc. 2000, 122, 11450. (22) Jurecka, P.; Sponer, J.; Cerny, J.; Hobza, P. Phys. Chem. Chem. Phys. 2006, 8, 1985. (23) Li, S.; Cooper, V. R.; Thonhauser, T.; Puzder, A.; Langreth, D. C. J. Phys. Chem. A 2008, 112, 9031. (24) Zimmerli, U.; Parrinello, M.; Koumoutsakos, P. J. Chem. Phys. 2004, 120, 2693. (25) Bludsky, O.; Rubes, M.; Soldan, P.; Nachtigall, P. J. Chem. Phys. 2008, 128, 114102. (26) Bludsky, O.; Rubes, M.; Soldan, P. Phys. ReV. B 2008, 77, 092103. (27) Rubes, M.; Bludsky, O. Phys. Chem. Chem. Phys. 2008, 10, 2611. (28) Rubes, M.; Bludsky, O.; Nachtigall, P. ChemPhysChem 2008, 9, 1702. (29) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (30) Soldan, P.; Hutson, J. M. J. Chem. Phys. 2000, 112, 4415. (31) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (32) Feyereisen, M.; Fitzgerald, G.; Komornicki, A. Chem. Phys. Lett. 1993, 208, 359. (33) Weigend, F. Phys. Chem. Chem. Phys. 2006, 8, 1057. (34) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297. (35) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (36) Gauss, J.; Stanton, J. F. Electron-correlated approaches for the calculation of NMR chemical shifts. AdV. Chem. Phys. 2002, 123, 355. (37) Blochl, P. E. Phys. ReV. B 1994, 50, 17953. (38) Ahlrichs, R.; Bar, M.; Haser, M.; Horn, H.; Kolmel, C. Chem. Phys. Lett. 1989, 162, 165. (39) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03; Gaussian, Inc., Wallingford, CT, 2004. (40) Molpro is a package of ab initio programs written by Amos, R. D.; Bernhardsson, A.; Berning, A.; Celani, P.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Hampel, C.; Hetzer, G.; Knowles, P. J.; Korona, T.; Lindh, R.; Lloyd, A. W.; McNicholas, S. J.; Manby, F. R.; Meyer, W.; Mura, M. E.; Nicklaß, A.; Palmieri, P.; Pitzer, R.; Rauhut, G.; Schu¨tz, M.; Schumann, U.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Werner, H.-J. (41) Kresse, G.; Hafner, J. Phys. ReV. B 1993, 48, 13115. (42) Suzuki, S.; Green, P. G.; Bumgarner, R. E.; Dasgupta, S.; Goddard, W. A.; Blake, G. A. Science 1992, 257, 942. (43) Gregory, J. K.; Clary, D. C. Mol. Phys. 1996, 88, 33. (44) Zhao, Y.; Tishchenko, O.; Truhlar, D. G. J. Phys. Chem. B 2005, 109, 19046. (45) Courty, A.; Mons, M.; Dimicoli, I.; Piuzzi, F.; Gaigeot, M.-P.; Brenner, V.; de Pujo, P.; Millie, P. J. Phys. Chem. A 1998, 102, 6590. (46) Cheng, B. M.; Grover, J. R.; Walters, E. A. Chem. Phys. Lett. 1995, 232, 364. (47) Jenness, G. R.; Jordan, K. D. J. Phys. Chem. C. Submitted for publication.
JP901410M
