Insights on the Role of Many-Body Polarization Effects in the Wetting

Subscriber access provided by READING UNIV

Article

Insights on the Role of Many-Body Polarization Effects in the Wetting of Graphitic Surfaces by Water Rahul Prasanna Misra, and Daniel Blankschtein J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.7b08891 • Publication Date (Web): 12 Nov 2017 Downloaded from http://pubs.acs.org on December 1, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry C is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Insights on the Role of Many-Body Polarization Effects in the Wetting of Graphitic Surfaces by Water Rahul Prasanna Misra and Daniel Blankschtein* Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA *Corresponding Author Daniel Blankschtein Herman P. Meissner (1929) Professor of Chemical Engineering Department of Chemical Engineering Massachusetts Institute of Technology Room 66-442B 77 Massachusetts Ave Cambridge MA 02139 USA email: [email protected] phone: 617.253.4594 fax: 617.253.4594

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 45

ABSTRACT It is well-known that atoms in a substrate placed in contact with a polar solvent like water experience a finite electric field from the solvent molecules. Nevertheless, the effect of this electric field on the wetting properties of the substrate remains unknown. In this study, by carrying our molecular dynamics (MD) simulations with force field parameters derived from ab initio simulations, we develop a theoretical framework to quantify the role of the polarization of graphene in the wetting of graphitic surfaces by water. Our study shows that a self-consistent modeling of the polarization of graphene yields a water contact angle on graphite which is remarkably different from the contact angle which results if the polarization energy is instead modeled implicitly using a Lennard-Jones potential, a typical approximation used in all previous MD simulation studies on the wetting of graphitic surfaces. Our findings reveal that polarization has a more pronounced effect on the interfacial entropy of water compared to dispersion. Consequently, polarization and dispersion effects contribute differently to the wetting of graphitic surfaces. Our study significantly advances our understanding of the water-graphene interface which is important for practical applications of graphene-based nanomaterials in osmotic power harvesting and seawater desalination.


Page 3 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


INTRODUCTION

Graphene, an allotrope of carbon and the two-dimensional form of bulk graphite, is known to possess extraordinary thermal, mechanical, and electronic properties. As a result of its atomic thickness and high mechanical strength, graphene has attracted widespread attention for use in potential membranes for seawater desalination,1,2 nanofluidic energy harvesting,3 nanoporebased DNA sequencing,4–6 and as a model system to explore the unique phase behavior of nanoconfined water.7,8 Although graphene has been a subject of intense research in the past decade, we still lack a thorough understanding of the mechanisms governing the interactions of water molecules with graphene. A promising way to characterize the nature and strength of these interactions is to develop a deeper understanding of how water wets graphene. In the past, the wetting transparency of graphene, that is, the extent to which graphene distorts the wetting properties of substrates placed below it, has elicited great debate, with reported claims of complete,9 partial,10,11 and negligible12 wetting transparency. In addition to a number of previous computational and experimental studies9–12 which investigated the wetting transparency of graphene, the wetting behavior of multilayer graphene and bulk graphite has also attracted significant attention in the scientific community. Past experimental measurements13,14 suggested that graphite is hydrophobic, with a water contact angle in the range of 84-86 ° . However, recent experimental data suggest that: (1) graphite is mildly hydrophilic, with a water contact angle in the range of 61 − 68° ,15–18 and (2) the higher contact angles measured in earlier experimental studies13,14 were due to hydrocarbon contamination on the graphite surface. In a recent study, Kozbial et al.19 studied the contact angle of water on freshly exfoliated graphite of varying qualities, and reported that the static contact angle of water on defect-free graphite is 65° ± 2°. It is important to stress that all previous studies on the transparency of graphene to van der Waals



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 45

interactions, as well on the wetting of multilayer graphene and graphite, have assumed that the water molecules interact with the carbon atoms in graphene solely through weak, pair-wise additive dispersion interactions. However, when the atoms of a substrate like graphite come in contact with a polar solvent like water, the permanent multipoles (charges, dipoles, quadrupoles, etc.) of the molecules in the polar solvent exert a finite electric field on the atoms of the substrate. This finite electric field can in turn distort the charge distribution in the substrate, resulting in induced charges whose magnitudes depend on the dielectric properties of the substrate. The resulting interaction energy between the permanent multipoles of the solvent molecules and the induced atomic multipoles in the substrate, referred to as the induction or polarization energy,20 is strictly pairwise nonadditive. Indeed, the electric field felt by an atom in the substrate is a vector whose magnitude and direction depend on the electric field exerted by all the solvent molecules and all the atoms in the substrate. Although the role of polarization effects on the structure and dynamic properties of water at a solid-liquid interface21–23 has been investigated in the past, to the best of our knowledge, the effect of many-body polarization effects on the contact angle of a polar solvent like water on any substrate has not yet been investigated. Note that ab initio molecular dynamics (AIMD) simulations, where the atomic motion is governed by intermolecular forces, determined from electronic structure calculations, can automatically capture the polarization effects described above, and thereby provide a valuable tool to study interfacial phenomena at the graphene-water interface.24 However, AIMD simulations are not suitable to calculate the contact angle of water on graphene by simulating water nanodroplets in contact with a graphene layer primarily because of: (i) significant line tension effects which can arise due to the inability of AIMD simulations to simulate systems 4 ACS Paragon Plus Environment

Page 5 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


containing more than a few hundred water molecules,25 (ii) the requirement of a much higher temperature (330-360 K) for AIMD simulations26 carried out using density functional theory with dispersion corrections (e.g., DFT+D) to model water in a liquid state, and (iii) the inability of AIMD simulations carried out using DFT with Grimme’s dispersion corrections27,28 to reproduce the experimental surface tension of bulk water at ambient conditions,29 which in turn can adversely affect the simulated contact angle.30 In a recent study, Nagata et al.29 calculated the surface tension of water at the water-air interface using AIMD simulations carried out using DFT with Grimme’s D2 and D3 dispersion corrections. Nagata et al.29 reported that the simulated values of the surface tension of water depend critically on the basis sets used, with simulations carried out using the double- valence polarized (DZVP) basis set at the BLYP/DZVP-D3 level of theory resulting in a dramatic overprediction of the surface tension of water (152 mJ/m2 instead of the experimental value of 72 mJ/m2). While the use of more complete basis sets like the triple- and quadruple- results in a remarkable improvement in the simulated surface tension of water, the surface tension of water predicted at the BLYP/TZV2P-D3 level of theory by Nagata et al.29, including the contribution due to finite size effects, still overpredicts the experimental value by ~35 mJ/m2. By using classical MD simulations, Leroy et al.30 reported that a difference of around 15 mJ/m2 in the simulated values of the surface tension of water using the TIP3P31 and TIP4P/200532 water models could, in turn, result in the simulated contact angles of water nanodroplets on graphene predicted by these water models differing by ~30°. Therefore, the inability of AIMD simulations to precisely reproduce the surface tension of water at the airwater interface at ambient pressure and temperature could potentially also affect the contact angles estimated from AIMD simulations of water nanodroplets in contact with a graphene sheet.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 45

We note here that an alternate methodology to estimate the contact angle of water on graphene is to first estimate the work of adhesion of water on graphene using a thermodynamic approach (e.g., the thermodynamic integration method30,33,34 or the free-energy perturbation method35), and then relate the simulated work of adhesion to the contact angle using the Young–Dupré equation36 given by: WSL = γ L (1 + cos θ )

(1)

where WSL is the work of adhesion of water (liquid) on graphene (solid), γ L is the surface tension of water, and θ is the contact angle of water on graphene. Note that the solid-liquid work of adhesion is defined as the reversible change in the Helmholtz free energy per unit area required to separate unit areas of solid and liquid from molecular contact to infinity (where they no longer interact). The work of adhesion, therefore, has energetic and entropic components. For water molecules interacting with graphene through a Lennard-Jones potential, Leroy et al.30 showed that the water-graphene work of adhesion is insensitive to the water model utilized to carry out the force field-based MD simulations, and does not depend on whether the water model used can correctly reproduce the surface tension of water. Therefore, according to the study by Leroy et al.,30 the experimental value of the surface tension of water can be used in conjunction with the obtained value of the work of adhesion from MD simulations to estimate the contact angle of water on graphene using Eq. (1). Following a similar approach, if the work of adhesion of water on graphene is estimated from AIMD simulations, then reliable estimates for the contact angle of water on graphene can be obtained using Eq. (1), even if the AIMD simulations are not able to precisely reproduce the surface tension of water. We note here that although the energetic component of the water-graphene work of adhesion can be readily obtained from AIMD simulations, the entropic contribution is difficult to estimate from AIMD simulations. Indeed, 6 ACS Paragon Plus Environment

Page 7 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


unlike force field-based MD simulations where the intermolecular interactions between the water molecules and graphene can be tuned, the intermolecular interactions in AIMD simulations are calculated from first-principles at every time step of the simulation, and therefore, it is not straightforward to apply thermodynamic approaches to estimate the work of adhesion from AIMD simulations, as done previously in several force field-based MD simulation studies.30,34,37,38 Moreover, estimating the work of adhesion by progressively moving the water film away from graphene, including calculating the interaction energy at different separation distances between the water film and graphene, using AIMD simulations would be practically unfeasible due to the very large computational cost for AIMD simulations. Therefore, an optimal methodology to investigate the wetting of graphitic surfaces is to first determine the force field parameters for MD simulations from highly accurate ab initio simulations which can then be used to carry out force field-based MD simulations. Force field-based MD simulations are still the method of choice to investigate the wetting of graphitic surfaces because (i) simulations on thousands of water molecules can be carried out efficiently for several nanoseconds of simulation time, thereby surpassing the length and time scale constraints in AIMD simulations, and (ii) some of the water models like the polarizable SWM4-NDP39 water model can predict water surface tension, and bulk water density, dielectric constant, viscosity, and diffusion coefficient values which are in good agreement with the experimental values at ambient conditions. In this study, by carrying out molecular dynamics (MD) simulations, using force field parameters derived from ab initio simulations, we develop the first theoretical framework to study the effect of the polarization of graphene on the wetting of graphitic surfaces by water. We model the induction energy contribution to the binding energy of a water molecule with graphene self-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 45

consistently using the classical Drude oscillator model proposed by Lamoureux and Roux,40 with Thole damping.41 We do this to answer the following long-standing, yet critical, question --- Is there a quantitative difference in the value of the simulated contact angle using an explicit and self-consistent treatment of polarization effects instead of using the conventional approach in previous MD simulation studies where the interaction energy between water molecules and graphene has been modeled using a pair-wise additive Lennard-Jones (LJ) potential? To answer this question, we first develop a theoretical framework to estimate the static dipole polarizability of the carbon atoms in graphene, and then derive an expression for the induction energy of the system consisting of water molecules and graphene. The derived expression for the induction energy is subsequently used to calculate the work of adhesion of water on graphite using the free-energy perturbation (FEP) method,35 which is then used to estimate the contact angle of water on graphite using the Young–Dupré equation.36 Our study, for the first time, points to a direct correlation between the electric field exerted by the water molecules on the carbon atoms in graphene and the orientation of the water molecules in the interfacial region. The electric field is, in turn, found to strongly depend on the static dipole polarizability of the carbon atoms in graphene. We also show that induction has a more pronounced effect on the interfacial entropy of water compared to dispersion, thus highlighting the necessity to model the induction contribution self-consistently in future studies of wetting-related interfacial phenomena of water on graphite or any graphene-coated substrate.

THEORETICAL METHODS

Ab initio Simulations. The binding energy of a water molecule with a graphene sheet was estimated by extrapolating the binding energy of water molecules interacting with very large polycyclic aromatic hydrocarbons (PAHs) like the circumcoronene (C54H18) and the dicircumcoronene (C96H24) molecules using the symmetry-adapted perturbation theory (SAPT).42 8 ACS Paragon Plus Environment

Page 9 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In this study, we used the simplest of the SAPT methods, SAPT0, where the monomers are treated at the Hartree-Fock level and electron correlation effects to dispersion interactions are estimated to second-order using perturbation theory. The binding energy of two interacting monomers using the SAPT0 method43 is given by: (10) (10) (20) (20) (2) (20 )  +  Eexch  +  Eind   (20)  E SAPT 0 =  Eelst + Eexch −ind + δ EHF  ind +  Edisp + Eexch − disp  disp elst exch

(2)

where E SAPT 0 is the binding energy which is decomposed into electrostatics, exchange, induction, and dispersion contributions. Additional discussion on each of the terms in Eq. (2) is included in Section S1 of the Supporting Information and also explained in detail in the study by ij Parker et al.43. Each of the terms in Eq. (2) is of the form E , where i represents the order of the

perturbation correction, and j represents the order for the intramonomer correlation in each of the monomers. Note that j = 0 in the SAPT0 method, indicating that intramonomer correlation effects are neglected. Parker et al.43 recently reported a new formulation for SAPT0, referred to (20) (20) as sSAPT0, by moderately modifying the SAPT0 Eexch −ind and Eexch − disp values which for certain

strongly interacting systems can correct for the strong overbinding at short range obtained using the SAPT0 method. When used in conjunction with the truncated double- ζ basis set - jun-ccpVDZ, Parker et al.43 reported obtaining very accurate binding energies using sSAPT0/jun-ccpVDZ for several databases of molecules like the S22 (e.g., water dimer) and NBC10 (e.g., benzene dimer), when compared with CCSD(T)/CBS data (widely considered to be the gold standard in quantum chemistry). Accordingly, Parker et al.43 recommended using sSAPT0/juncc-pVDZ as the bronze standard for SAPT in order to calculate intermolecular interaction energies. The ab initio simulations used in this study for estimating the binding energies of water molecules with PAHs were carried out using the sSAPT043 method with a jun-cc-pVDZ basis set 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 45

as implemented in the psi444 open source program. In our study, we used the sSAPT0 method because of the significantly lower computational cost of sSAPT0 compared to other SAPT techniques, while still maintaining a high level of accuracy which allowed us to calculate the binding energy of a water molecule with larger PAHs like the dicircumcoronene molecule (C96 H 24 ) . This allowed us to study the convergence of various contributions to the binding energy as the size of the PAHs was increased to approach the limit of an infinitely-periodic graphene lattice. Molecular Dynamics Simulations. The molecular dynamics (MD) simulations reported in our study were carried out using the LAMMPS software package45 in the canonical NVT (constant number of particles, volume, and temperature) ensemble. In the MD simulations carried out to estimate the orientation of interfacial water molecules and for calculating the work of adhesion of water on graphite, we considered a water film consisting of 1000 water molecules interacting with graphite. Graphite was modeled using five graphene layers with the carbon atoms in the top three graphene layers allowed to vibrate, while the carbon atoms in the bottom two graphene layers were kept fixed to support the top three graphene layers. We note here that a similar approach to keep the bottom layers fixed when a substrate is brought in contact with a graphene layer was used in an earlier study.46 However, in order to reduce computational cost, we neglected all LJ interactions between the water molecules and the fixed carbon atoms in the fourth and fifth graphene layers. Note that the carbon atoms in the third graphene layer contribute minimally to both the dispersion and polarization energies, which allowed us to model the interaction of water molecules with bulk 3D graphite using solely the top three graphene layers which interact with the water molecules (see Section S5 in the Supporting Information document for additional details). 10 ACS Paragon Plus Environment

Page 11 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Our periodic simulation box contained 1000 water molecules (modeled using the polarizable SWM4-NDP water model39) interacting with graphite, and the dimensions in the x, y, and z directions were 27.02 Å , 25.53 Å , and 200 Å , respectively. The dimension in the z direction was intentionally kept large to prevent the interaction of periodic images. We used FFTool47 and Packmol48 utilities to prepare the initial configurations of our system. The carbon atoms in the top three graphene layers were modeled using the Dreiding force field,49 which can reproduce the phonon dispersion relation of graphene with reasonable accuracy.50 In the MD simulations where we considered the polarization of graphene from the electric field exerted by the water molecules, the induced dipoles on the carbon atoms were modeled using the classical Drude oscillator model40 with Thole damping.41 In the classical Drude oscillator model, an induced dipole is modeled as a pair of Drude particle (DP) and Drude core (DC) possessing point charges, with the DP attached to the DC by a harmonic spring, and an extended dynamical approach is followed where the DPs are maintained at a significantly lower temperature (1K), compared to the DCs, which are maintained at close to room temperature (300 K). In our simulations, the DPs attached to the carbon atoms were assigned a small mass,

mD = 0.4 g mol−1 ,

the

carbon

atom

which

is

the

DC

was

assigned

a

mass,

mC = 12.011 − 0.4 = 11.611 g mol−1 , and two separate Nosé-Hoover thermostats51 were used to maintain the temperature of the DPs and DCs. The temperature of the DCs, i.e., of the carbon atoms in the top three graphene layers, was maintained at 300 K using a Nosé-Hoover thermostat with a time constant of 0.1 ps. The temperature of the DPs attached to the carbon atoms in graphene was maintained at 1 K using a separate Nosé-Hoover thermostat with a time constant of 5 fs. All the MD simulations were carried out using a time step of 1 fs. A direct-space cutoff of 15 Å was used for both the LJ and the Coulombic interactions and the long-range part of the 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 45

Coulombic interactions was calculated using the PPPM method52 with a tolerance of 1 × 10 −4 in the energy. We note here that previous studies30,53,54 have shown that the simulated contact angles obtained from MD simulations can be sensitive to the cutoff length used to truncate the LJ interactions. For this reason, we have utilized a large LJ cutoff length of 15 Å to carry out the work of adhesion simulations reported in this study. To test the effect of the LJ cutoff length on our results, we carried out an additional MD simulation of a film of water interacting with graphite, utilizing an MD simulation setup which is identical to the one that we used to carry out the work of adhesion simulations. Upon using a LJ cutoff length of 20 Å, the deduced value of the energetic component of the work of adhesion differed by less than 3% from the result obtained using a 15 Å cutoff. This indicates that our predictions using a large LJ cutoff length of 15 Å are reasonably accurate, and that use of a higher LJ cutoff length will not significantly change any of the results presented in this study.

RESULTS AND DISCUSSION

Water-Graphene Interaction. First, we follow a methodology similar to the one introduced by Jenness et al.55,56 to estimate the binding energy of a water molecule interacting with graphene. With the goal of estimating the induction energy contribution to the binding energy, we carried out ab initio simulations using the sSAPT043 method with a jun-cc-pVDZ basis set as implemented in the psi444 open source program. The sSAPT0 method decomposes the binding energy into four physically meaningful contributions: electrostatics, exchange, dispersion, and induction. We estimated the binding energy of a water molecule interacting with very large polycyclic aromatic hydrocarbons (PAHs) like the circumcoronene (C54H18) and the dicircumcoronene (C96H24) molecules using the sSAPT0 method and show that the binding energy of a water molecule with large PAHs like the circumcoronene molecule can be


Page 13 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


conveniently extrapolated to obtain the binding energy of a water molecule with a periodic graphene layer (see Section S1 in the Supporting Information document for additional details). In order to test the accuracy of our predictions using the sSAPT0 method, we first compared our results for the binding energy of a water molecule with benzene against highly accurate coupled cluster with single, double, and perturbative triple excitations (∆CCSD(T)) , and with quantum Monte Carlo (DMC) results reported by Ma et al.57. The binding energy curve obtained using the sSAPT0 method is in reasonable agreement with the ∆CCSD(T) and DMC results. For example, compared to the ∆CCSD(T) and DMC results, the error in the magnitude of the binding energy minimum obtained using the sSAPT0 method is less than 4% (see Figure S1 in the Supporting Information document). Next, we compared our predictions for the binding energy of a water molecule with a periodic graphene layer obtained using the sSAPT0 method with DFT simulations carried out by Al-Hamdani et al.58 using the dispersion corrected exchangecorrelation functional: PBE+D3,28,59 and dispersion inclusive exchange correlational functionals: vdW-DF60,61 and vdw-DF262. We found that in comparison to the results from DFT simulations, our predictions for the binding energy using the sSAPT0 method are in better agreement with benchmark results by Ma et al.63 obtained using the highly accurate random-phase approximation (RPA) method (see Figure S5 in the Supporting Information document). Furthermore, we have also carried out additional test simulations to study the effect of the lateral position of a water molecule placed above a circumcoronene molecule on the binding energy, and found that the various contributions to the binding energy do not vary appreciably based on the lateral position of the water molecule (difference in the total binding energy < 0.03 kcal/mol), a finding that is consistent with a previous study by Tocci et al.24. The interested reader is referred to Section S1



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 45

in the Supporting Information document for additional details on these test simulations and comparison to benchmark results. We note here that by using the sSAPT0 method, the electrostatic energy of interaction of a water molecule with a circumcoronene molecule can be decomposed into a short-range charge penetration contribution (arising from the overlap of the charge distributions of interacting molecules20) and a long-range permanent multipole moments contribution, using the Gaussian Distributed Multipole Analysis (GDMA) technique.64 Following previous studies,53,65 it can be shown that only the short-range charge penetration contribution to the electrostatic energy contributes towards the binding energy of a water molecule with a periodic graphene layer. Because the short-range charge penetration and exchange contributions to the binding energy are pair-wise additive,20 while the dispersion contribution is also approximately pair-wise additive,20 we model these three contributions to the binding energy together using a pair-wise additive Lennard-Jones (LJ) potential, U LJ

 σ co 12  σ co 6  = 4ε co   −   , where U LJ is the LJ interaction  r    r 

potential of a water molecule with a carbon atom in graphene, r is the separation distance between the oxygen atom of a water molecule and a carbon atom in graphene, and ε co and σ co are the LJ parameters determined by fitting to the sum of the charge-penetration, exchange, and dispersion contributions to the binding energy obtained from ab initio simulations carried out using the sSAPT0 method. An important consideration in fitting the LJ parameters, ε co and σ co , to model the charge penetration, exchange, and dispersion contributions to the binding energy is to choose the proper orientation of the water molecule to fit the LJ parameters. Note that the LJ interaction potential is spherically-symmetric, which makes it difficult to model the dispersion energy of water 14 ACS Paragon Plus Environment

Page 15 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molecules in various orientations using a single value for the parameters, ε co and σ co . The binding energy plots in Figure 1 (a) clearly indicate that for a single water molecule interacting with a circumcoronene molecule, the 2 leg and 1 leg orientations are lower in energy, i.e., more stable relative to the tangential and inverted 1 leg orientations. In Figure 1 (b), we show the fitted LJ potential for a tangentially-oriented water molecule against reference sSAPT0/jun-cc-pVDZ simulations. We have also fitted the LJ parameters using reference sSAPT0/jun-cc-pVDZ simulations obtained for the 1 leg orientation. Next, we carried out MD simulations of a water film consisting of 1000 water molecules interacting with graphite with the LJ potential fitted to the tangential and 1 leg orientations, respectively (see Supporting Information, Section S1 for additional details on the MD simulation setup). As shown in Figure 1 (c) and Figure 1 (d), we found that irrespective of the orientation of the water molecule used to fit the LJ parameters against the reference sSAPT0/jun-cc-pVDZ simulations, the dominant orientation of the water molecules in the first hydration layer (defined as the region from the surface of the 1st graphene layer until the location of the 1st minimum in the density profile of water molecules) is always tangential. Indeed, we found that water molecules in the first hydration layer prefer to remain in the tangential and inverted 1 leg orientations in approximately a 2:1 ratio, and that this is largely independent of the orientation of the water molecule chosen to fit the LJ parameters. This follows because the water molecules near the graphite interface seek to maximize the number of hydrogen-bonded neighbors by remaining in plane (tangential orientation), while the 2 leg and 1 leg orientations result in the loss of hydrogen-bonds for those –OH bonds which point towards the graphite surface. Hence, we show that the dominant orientation of the interfacial water molecules in a water film placed over graphite can be different from the most stable water orientations which result when a single isolated water molecule interacts with graphene. Note



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 45

that in the MD simulations discussed above, we have not yet considered the polarization of the carbon atoms in graphene. However, in Figure S12 (b) in the Supporting Information document, we show that the dominant orientation of the water molecules in the interfacial layer remains tangential even when we consider the polarization of graphene in the presence of the water molecules. Moreover, using surface-sensitive sum frequency generation (SFG) spectroscopy, Singla et al.66 recently reported that the orientation of the water molecules next to graphene is planar (tangential), an experimental result which is in good agreement with the results from the MD simulations presented here. Having determined that the dominant orientation of water molecules in the interfacial layer is tangential, the LJ parameters, ε co and σ co , were fitted to the sum of charge-penetration, exchange, and dispersion components of the binding energy from sSAPT0/jun-cc-pVDZ simulations of the interaction of a tangentially-oriented water molecule with circumcoronene. The values of the fitted LJ parameters are: ε co = 0.1072 kcal/mol and σ co = 3.3215 Å . Because the induction energy component of the binding energy is strictly pair-wise non-additive,20 we model this component of the binding energy self-consistently using the classical Drude oscillator model.40 In previous MD simulation studies on the wetting of graphitic surfaces, the total binding energy, including the induction energy, was modeled using a Lennard Jones potential. It is noteworthy that it is possible to fit the LJ parameters to model the induction and dispersion energies together using a LJ potential, because the leading-order terms in both the induction (permanent dipole – induced dipole interaction) and the dispersion (instantaneous dipole – induced dipole interaction) energies scale as 1/ r 6 , where r is the separation distance between a water molecule and a carbon atom in graphene. To compare the contact angle of water on graphite obtained using an implicit modeling of polarization effects using a LJ potential (an 16 ACS Paragon Plus Environment

Page 17 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


unsubstantiated approximation used in previous MD simulation studies) with the case where the induction energy is modeled self-consistently using the classical Drude oscillator model, we derive a second set of LJ parameters, ε co = 0.1448 kcal/mol, and σ co = 3.2148 Å, obtained by fitting the LJ parameters to reproduce the total binding energy of a tangentially-oriented water molecule with circumcoronene obtained from sSAPT0/jun-cc-pVDZ simulations (see Figure 1(b)).

Figure 1. (a) Binding energy plots of a water molecule interacting with a circumcoronene molecule obtained using sSAPT0/jun-cc-pVDZ simulations for different water orientations. The water-circumcoronene distance is defined as the vertical distance between the oxygen atom in the water molecule and the center of the circumcoronene molecule. Side view of the various orientations shown here correspond to (i) 1 leg, (ii) 2 leg, (iii) tangential, and (iv) inverted 1 leg orientations, respectively. The inset shows the top view of a tangentially-oriented water molecule interacting with a circumcoronene molecule. (b) Fitted LJ potential which can reproduce the interaction energy corresponding to the sum of the charge-penetration, exchange, and dispersion contributions to the binding energy obtained from sSAPT0/jun-cc-pVDZ simulations (denoted by green diamonds), and for the case where the LJ parameters are fitted to reproduce the total binding energy, including the induction energy from sSAPT0/jun-cc-pVDZ simulations (denoted by red circles). The binding energy curves shown here are for a tangentially-oriented water 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 45

molecule. (c) and (d) Probability distribution of cos(φ ) where φ is the angle which the plane of a water molecule makes with the plane of graphite. The LJ parameters, ε co and σ co , for the MD simulations in (c) and (d) were fitted to reproduce the sum of the charge-penetration, exchange, and dispersion contributions to the binding energy from sSAPT0/jun-cc-pVDZ simulations of a water molecule interacting with a circumcoronene molecule in the 1 leg and tangential orientations, respectively.

Modeling the Polarization of Graphene. For the purpose of modeling the polarization of the carbon atoms in graphene in the presence of water molecules, we first develop a theoretical framework to estimate the static dipole polarizability of the carbon atoms in graphene. The polarizability of multilayer graphene films under an applied electric field has been studied by Yu et al.67 using Density Functional Theory (DFT), with a plane wave basis set within the local density approximation (LDA). The DFT-LDA calculations of Yu et al.67 based on a Berry phase approach predict a linear increase in the polarizability of graphene films in the out-of-plane direction (α C⊥ ) with the number of graphene layers. Further, based on the study of the dielectric and optical properties of graphite, Phillips68 recommended an anisotropy ratio for the components of the polarizability tensor in the in-plane (α C|| ) and out-of-plane (α C⊥ ) directions, with

αC

||

αC

= 3.5. We note here that similar to graphene films, carbon nanotubes (CNTs) also

⊥

show an anisotropy in the polarizability along the axial and transverse directions.69,70 In order to model the anisotropy in the components of the polarizability tensor of graphene in the in-plane and out-of-plane directions which is a challenging endeavor in MD simulations, we utilized the classical Drude oscillator model40 in conjunction with Thole’s damped interaction model,41 resulting in two parameters for our model: the static dipole polarizability of the carbon atoms in graphene ( α C ) and the Thole damping parameter ( aC ).41 p q For the interaction between two dipoles, µ and µ , the components of the dipole field tensor


Page 19 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


following the notation by Stone,20 is given by: 2 1  3Rα Rβ − R δαβ =  4πε 0  R5

pq

Tαβ

  

(3)

i i i where {α , β } ∈ { x, y , z} , {Rα , Rβ } ∈ {x p − x q , y p − y q , z p − z q } , x , y , z are the positions of the

i th

induced

dipole

in

the

x,

y,

and

z

directions,

respectively,

1 if α = β  R = ( x p − x q ) 2 + ( y p − y q ) 2 + ( z p − z q ) 2 , and δαβ =   . In this work, we have used 0 if α ≠ β  the exponential form of the Thole dipole field tensor (this corresponds to the 1st entry in Table 1 in Thole’s study41) which is given by:

Tαβpq =

 3uα uβ  5 α pα q  u

1 4πε 0

where ui =

Ri

(α α )

1/6

p

q

, u=

  a 3u 3 a 2 u 2  δ  + + au + 1 e − au  − αβ3 1 −  2    6  u

R

(α α )

1/6

p

  a 2u 2   + au + 1 e − au   (4) 1 −      2 

, αi stands for the dipole polarizability of the i th atom in

q

units of Å3 , and a is the dimensionless Thole damping parameter. Note that higher values of a indicates lesser damping, and one recovers the original dipole field tensor in the limit of a → ∞ , while a =0 results in a scenario where the induced dipoles do not interact with each other, i.e.,

Tαβpq = 0. Next, we consider a unit cell containing two carbon atoms, C1 and C2 (shown in Figure 2) in a periodic graphene lattice. Let E represent the electric field felt by the carbon atoms, C1 and C2 2 in the unit cell, and µ1ind and µind represent the induced dipole moments at the positions of the

carbon atoms, C1 and C2 . Because the induced dipoles interact with each other in the periodic 2 graphene lattice, the induced dipole moments, µ1ind and µind , are given by20:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1 ind

µ

µ

2 ind

where r1 =

Page 20 of 45

  ∞ ∞ ∞ ∞ 11 1 12 2   = α C E + ∑ ∑ T ( r1 , r1 + ia1 + ja 2 ) .µ ind + ∑ ∑ T ( r1 , r2 + ia1 + ja 2 ) .µ ind   i =−∞ j =−∞ i =−∞ j =−∞   j =i ≠ 0

(5)

  ∞ ∞ ∞ ∞ 21 1 22 2   = α C E + ∑ ∑ T ( r2 ,r1 + ia1 + ja 2 ) .µ ind + ∑ ∑ T ( r2 , r2 + ia1 + ja 2 ) .µ ind   i =−∞ j =−∞ i =−∞ j =−∞   j =i ≠ 0

(6)

1 2 l 3 ( a1 + a 2 ) , r2 = ( a1 + a 2 ) , a1 = lxˆ , a 2 = xˆ + lyˆ , αC is the dipole polarizability 3 3 2 2

tensor whose diagonal elements are equal to αC and all the off-diagonal elements are equal to zero, l is the lattice constant for graphene equal to 2.46 Å , a1 and a2 are the lattice vectors for the unit cell, r1 and r2 are the position vectors of the carbon atoms in the unit cell described in j terms of the lattice vectors, Tij (ri ,r j ) . µ ind represents the electric field at the i th induced dipole

th located at position ri from the j induced dipole located at position r j , and the individual

components of the dipole field tensor, Tij (ri ,r j ) , are calculated using Eq. (4). Note that from the definition of the dipole field tensor, Tij , in Eq. (4), T12 = T21 , and T11 = T22 (from symmetry), 2 which implies that, µ1ind = µ ind from Eqs. (5) and (6). The net induced dipole moment of the unit

cell is then given by: −1

cell µind = 2µ1ind

  ∞ ∞ ∞ ∞ −1 11 12  = 2 α C − ∑ ∑ T ( r1 , r1 + ia1 + ja 2 ) − ∑ ∑ T ( r1 , r2 + ia1 + ja 2 )  E = α cell E (7)   i =−∞ j =−∞ i =−∞ j =−∞ j =i ≠ 0  

Although the summations in Eq. (7) extend until infinity, we obtained well-converged solutions

( error in α

cell

αβ

< 10 −3 Å 3 ) by considering 2000 × 2000 supercells in the x and y directions,

respectively. We stress here that due to the symmetry of the periodic graphene lattice, from our model, the off-diagonal elements of α cell are equal to zero, and α xxcell = α yycell ≠ α zzcell . Note that the 20 ACS Paragon Plus Environment

Page 21 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


magnitudes of α xxcell and α zzcell depend solely on the two parameters, αC and aC , which we determine uniquely from Eq. (7) by requiring that α zzcell = 0.867 Å3 ( Yu et al.67 have reported,

α

cell zz

cell α xxcell α yy = 0.867 Å for single-layer graphene using DFT-LDA simulations), and cell = cell = 3.5 α zz α zz 3

(Phillips68 recommended this ratio based on dielectric and optical studies of graphite as discussed earlier which has also been used in some previous studies71,72). By solving Eq. (7) numerically using these two conditions, we obtain the polarizability parameters, α C = 1.139 Å3 and a C = 1.507 , which can reproduce the polarizability tensor of a periodic graphene lattice (see Section S2 in the Supporting Information for details on additional calculations carried out for multilayer graphene films).

Figure 2: Schematic of a periodic graphene lattice with C1 and C2 carbon atoms marked in a unit cell used to estimate αc and ac . After estimating the polarizability parameters, α C and aC , which can correctly reproduce the components of the polarizability tensor of graphene, we next calculate the induction energy of a system consisting of polarizable water molecules interacting with polarizable carbon atoms in graphene. For reference, we use the 4-site SWM4-NDP39 water model which has an additional massless site, M, along the angle bisector of the H-O-H angle. The hydrogen atoms and the M



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 45

site have partial atomic charges which can adequately represent the gas-phase dipole and quadrupole moments of the water molecules, while the oxygen atom is polarizable with a dipole polarizability, α w . Following the approach20 that the induction energy resulting from the polarization of any atom is given by two contributions: (i) the interaction between the induced dipole (resulting from the polarization of the atom) with the permanent and induced dipoles of all the other atoms in the system, and (ii) the self energy of the induced dipole which is equal to the work required to distort the charge distribution in the atom to create the induced dipole, the total polarization energy of a system consisting of polarizable water molecules interacting with polarizable carbon atoms in graphene is given by:

U pol = −

1 1 % % %j i% j j µ iind . Tii . µ ind − ∑∑ µ ind . T jj . µ ind − ∑∑ µ iind . Tij . µ ind ∑∑ 2 i i% ≠i 2 j %j ≠ j i j

−∑∑ µ i

where E j =

∑

j

k = H1 , H 2 , M

i ind

. E − ∑∑ µ j

j

%j ≠ j

Ek , where E k =

j ind

1 1 . E + ∑ µ ind . α C −1 . µ ind + ∑ µ ind . αW −1. µ ind i i j j 2 i 2 j

(8)

%j

qk R is the electric field exerted by the partial atomic 4πε 0 R 3

k charge, q , in the hydrogen atoms and in the M site of the water molecules, and the indices i and

j are used to sum over all the induced dipoles in graphene and in the water molecules, respectively. We note here that the components of the dipole field tensor, T ii , describing the %

interactions between the induced dipoles in graphene, are given by the exponential form of the Thole dipole field tensor (see Eq. (4)), while the other dipole field tensors in Eq. (8) remain unchanged (given by Eq. (3)). Note that the first term in Eq. (8) reflects the interactions between the induced dipoles in graphene, the second term reflects the interactions between the induced dipoles in the water molecules, the third term reflects the interactions between the induced dipoles in graphene and the induced dipoles in the water molecules, the fourth term reflects the 22 ACS Paragon Plus Environment

Page 23 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


interactions between the induced dipoles in graphene and the partial atomic charges in the hydrogen atoms and the M site of the water molecules, the fifth term reflects the interactions between the induced dipoles in the water molecules and the partial atomic charges in the hydrogen atoms and the M site of the water molecules, the sixth term reflects the self energy of the induced dipoles in graphene, and the seventh term reflects the self energy of the induced dipoles in water molecules. Note that a factor of one half is included in the first and the second terms to correct for double counting. The values of the induced dipole moments are obtained by minimizing the energy, U pol . j Specifically, evaluating the gradient of U pol with respect to µiind and µind , and then setting these

to zero, we obtain expressions for the induced dipole moments in graphene ( µ iind ) and in the water molecules which are given by:

 1 1 % % j i% i% µ iind = α C  ∑ E j + ∑ Tij . µ ind + ∑ µ ind + ∑ Tii . µ ind . Tii  2 i% ≠i 2 i% ≠ i j  j 

(9)

  1 1 % % %j %j % j i = αW  ∑ E j + ∑ T jj . µind + ∑ µind µ ind . T jj + ∑ µ ind . Tij  2 %j ≠ j 2 %j ≠ j i  %j ≠ j 

(10)

( )

% T

Because the dipole field tensor is symmetric by definition,20 i.e., T jj

= T jj , it follows from the %

i j j rules of linear algebra that Tii . µiind = µind .Tii and T jj . µind = µind .T jj , such that Eqs. (9) and (10) can be rewritten as follows:

%

%

%

%

%

%

%

%

 % %  j µ iind = α C  ∑ E j + ∑ Tij . µ ind + ∑ Tii . µ iind  j i% ≠ i  j 

(11)

 % % %j j i  µ ind = αW  ∑ E j + ∑ T jj . µ ind + ∑ Tij . µ ind  %j ≠ j i  %j ≠ j 

(12)

j Substituting the expressions for µiind and µind from Eqs. (11) and (12), respectively, in Eq. (8), we

obtain a simplified expression for the polarization energy, U pol , which is given by:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

U pol = −

1 1 % . E j − ∑∑ µ ind . Ej µ ind ∑∑ i j 2 i j 2 j %j ≠ j

Page 24 of 45

(13)

Note that the polarization energy, U pol , in Eq. (13) consists of two contributions. The first one represents the polarization of the carbon atoms in graphene resulting from the electric field exerted by the water molecules, while the second one arises as a result of the polarization of the water molecules. In Figure 3 (a), we compare the induction energy (U pol ) calculated using the polarizability parameters, α C = 1.139 Å3 and a C = 1.507, using the classical Drude oscillator model, to the induction energy component of the binding energy of a tangentially-oriented water molecule with circumcoronene obtained using sSAPT0/jun-cc-pVDZ simulations. The interested reader is referred to Section S3 in the Supporting Information for additional details on the classical Drude oscillator model.

Figure 3. (a) Prediction of the induction energy of a tangentially-oriented water molecule interacting with circumcoronene obtained from sSAPT0/jun-cc-pVDZ simulations is compared to the induction energy calculated using the classical Drude oscillator model. (b) The induction energy of Na+ and Cl- ions interacting with circumcoronene obtained from SAPT0/jun-cc-pVDZ simulations is compared to the induction energy calculated using the classical Drude oscillator model.

Figure 3 (a) shows that the induction energy calculated using the classical Drude oscillator model underpredicts the induction energy obtained using sSAPT0/jun-cc-pVDZ simulations. This discrepancy may be attributed to charge penetration resulting from the overlap of the charge 24 ACS Paragon Plus Environment

Page 25 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


distributions of the water molecule and the carbon atoms in circumcoronene. Indeed, the possibility of a potential mismatch between the induction energy estimated using an ab initio perturbation theory, like SAPT, and any classical multipole-based approach due to short-range penetration effects around the region of the van der Waals minimum, has been reported previously by Jeziorski et al.73 To further support our hypothesis that the mismatch in the induction energy may be due to charge penetration effects, we carried out additional simulations using the SAPT0/jun-cc-pVDZ method to estimate the binding energies of the salt ions, Na+ and Cl-, interacting with circumcoronene. In this case, we expect that the magnitude of the induction energy should be larger than the one obtained for the water-circumcoronene interaction due to the high electric fields exerted by the salt ions.74 Our choice of studying the interactions of salt ions with circumcoronene is motivated by the fact that we know beforehand that polarization effects would be important for this system, thereby giving us another opportunity to test the predictions of the classical Drude oscillator model against reference ab initio data. We observed that in the case of Na+ and Cl- ions interacting with circumcoronene, the charge penetration contribution to the electrostatic energy rapidly decays to zero for separation distances greater than 3 Å , for both the Na+ and Cl- ions (see Figure S9 in the Supporting Information). Figure 3 (b) shows that for separation distances greater than 3 Å , where charge penetration effects are negligible, the induction energy for Na+ and Cl- ions interacting with circumcoronene, calculated using the classical Drude oscillator model (Eq. (13)), with polarizability parameters,

α C = 1.139 Å3 and aC = 1.507, is in excellent agreement with the result from ab initio simulations (see Supporting Information, Section S4 for additional details on the ab initio simulations carried out to estimate the binding energies of the Na+ and Cl- ions interacting with circumcoronene). It is worth mentioning here that, at a classical level of theory, the induction 25 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 45

energy of Na+ and Cl- ions interacting with circumcoronene should be identical. This is because the classical induction energy depends only on the magnitude of the charge of a salt ion and not on its sign. Indeed, we observe that the induction energy of Na+ and Cl- ions with circumcoronene obtained from ab initio simulations is identical for separation distances greater than 3 Å from the cirumcoronene molecule where quantum chemical charge penetration effects become negligible. On the other hand, when a water molecule interacts with circumcoronene, the charge penetration contribution to the electrostatic energy decays slowly and is not negligible near the minimum of the binding energy (see Figure S2 (d) in the Supporting Information), thereby supporting our hypothesis that charge penetration effects may be responsible for the mismatch in the induction energies calculated using the classical Drude oscillator model and from sSAPT0/jun-cc-pVDZ simulations. Since there exists no classical theory to model the induction energy resulting from quantum chemical charge penetration effects, we fitted the parameters, αC and aC , to the induction energy obtained from sSAPT0/jun-cc-pVDZ simulations. As shown in Figure 3 (a), the induction energy calculated with the fitted parameters, α C = 1.80 Å3 and a C = 0.50 , is in excellent agreement with the sSAPT0 data for all separation distances greater than 3 Å . Note that when carrying out MD simulations of water droplets on graphite, we expect that no water molecules should approach the top-most graphene layer closer than 3 Å due to the strong shortrange exchange repulsion from the LJ potential. Contact Angle of Water on Graphite. After determining the polarizability parameters, α C and

aC , we next investigate the role of polarization effects on the contact angle of water on graphite. For this purpose, we first estimate the work of adhesion of water on graphite ( WSL ), which is 26 ACS Paragon Plus Environment

Page 27 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


subsequently used to calculate the contact angle ( θ ) of water on graphite using the Young– Dupré equation36 (see Eq. (1)). For estimating WSL , we have used the free-energy perturbation (FEP) method of Zwanzig35 as follows:

 − (U1 − U 0 )  WSL = F1 − F0 = −k BT ln exp   k BT  

(14) 0

where F denotes the Helmholtz free energy of the system per unit area, U denotes the total interaction energy between the water molecules and graphite per unit area, subscript 0 is used to denote the initial state of the system, where the water molecules fully interact with graphite, subscript 1 is used to denote the final state of the system, where the water molecules no longer

 − (U1 − U 0 )  interact with graphite, exp   k BT  

is evaluated by sampling over the trajectory in state 0

0, k B is the Boltzmann constant, and T is the absolute temperature of the system. Note that U in Eq. (14) contains contributions from both the dispersion energy (U LJ ) and the induction energy

(U pol ). Because WSL is the reversible change in the Helmholtz free energy per unit area required to separate unit areas of graphite and water from molecular contact to infinity (where they no longer interact), we calculated WSL by devising a reversible path where we first tune αC from its actual value to zero, while keeping the parameter ε co in U LJ constant. We refer to this inductionrelated contribution to WSL as WSLind . Next, we tune the parameter ε co in U LJ , while maintaining

αC = 0 , from its actual value to zero. We refer to this dispersion-related contribution to WSL as disp WSLdisp . This two-step method allows us to conveniently calculate WSL = WSind L + WSL .



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 45

To estimate WSL using the methodology described above, we carried out MD simulations of a water film comprised of 1000 water molecules interacting with graphite. The MD simulations were carried out in the canonical NVT (constant number of particles, volume, and temperature) ensemble using the LAMMPS45 MD package. The water molecules were modeled using the polarizable SWM4-NDP water model,39 and the polarization of the carbon atoms in graphite was modeled using the classical Drude oscillator model,40 as implemented in the LAMMPS MD package by Dequidt et al.75 Note that we also carried out additional simulations using the nonpolarizable SPC/E76 water model and found that WSL is not sensitive to the water model used. We also found that the 2nd term in the expression for U pol in Eq. (13), which corresponds to the polarization of the water molecules, does not contribute to WSL due to energy-entropy compensation77 (see Supporting Information, Section S5 for a complete description of the computational setup used to carry out the MD simulations to calculate WSL ). In order to comprehensively study the roles of polarization and dispersion effects on the contact angle of water on graphite, it is sufficient to study four different cases: Case I, where the polarization of graphite is neglected and the water molecules interact with graphite solely through a LJ potential, with LJ parameters, ε co = 0.1072 kcal/mol and σ CO = 3.3215 Å , obtained by requiring that the LJ potential ( U LJ ) is able to model the pair-wise additive components of the binding energy of a water molecule with graphene, which include the short-range charge penetration, exchange, and dispersion energies, estimated from sSAPT0/jun-cc-pVDZ simulations; Case II, where the LJ potential is the same as in Case I, but we also consider the polarization of graphite which is modeled using the polarizability parameters, αC = 1.139 Å3 and a C = 1.507, which can reproduce the polarizability tensor of graphene; Case III, where the LJ 28 ACS Paragon Plus Environment

Page 29 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


potential is the same as in Cases I and II, and the polarizability parameters, αC = 1.80 Å3 and a C = 0.5, are obtained by fitting to the induction energy obtained using sSAPT0/jun-cc-pVDZ simulations, thereby also modeling the charge penetration contribution to the induction energy; and Case IV, where the LJ parameters, ε co = 0.1448 kcal/mol and σ co = 3.2148 Å , are obtained by fitting to the total binding energy of a water molecule interacting with graphene obtained using sSAPT0/jun-cc-pVDZ simulations. Note that for a single water molecule interacting with graphene, the force field parameters in Cases III and IV will result in the same binding energy. However, while Case III represents the correct methodology to self-consistently model the induction energy which is strictly pair-wise non-additive, Case IV represents just an approximate, implicit method to do so where the induction and dispersion energies are modeled together using a pair-wise additive LJ potential. Note that if Cases III and IV result in the same contact angle, then, this would indicate that the induction energy can be modeled using a LJ potential, an unsubstantiated assumption made in all previous MD simulations studies on the wetting of graphitic surfaces by water. In Table 1, we report the calculated WSL values at T=300 K

for

the

four

cases

discussed

above,

where

WSL = ∆U SL − T ∆SSL = WSLdisp + WSLind ,

disp disp ind ind disp ind disp ind WSLdisp = ∆U SL − T ∆SSL − T ∆S SL + ∆U SL + ∆SSL , WSLind = ∆U SL , ∆U SL = ∆U SL , ∆S SL = ∆SSL ,

∆U SL and T ∆S SL are the energetic and entropic contributions to WSL . By comparing Cases I and II in Table 1, we find that the contact angle of water on graphite decreases by a negligible amount (around 3°) when the polarization of graphite is considered and modeled using the polarizability parameters, α C and aC , which can reproduce the polarizability tensor of graphene. Note that the contact angle, θ , reported in the last column of Table 1 was estimated from W SL using Eq. (1) with γ L = 72 mJ/m 2 at T=300 K.78 For Case III, where α C 29 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 45

and aC are fitted to the induction energy obtained using sSAPT0/jun-cc-pVDZ simulations, the obtained value of the contact angle is 9° smaller than the one obtained in Case I. In contrast, we find that the contact angle calculated in Case IV is 32 ° smaller than the one calculated in Case I. Therefore, for the first time, our results show that the contact angle calculated using an implicit treatment of polarization effects (Case IV), is remarkably different (23 degrees) from the contact angle calculated using an explicit treatment of polarization effects (Case III). Table 1. Breakdown of WSL into induction and dispersion contributions, which are further decomposed into energetic and entropic contributions for different values of the interaction parameters.a

Case

Interaction Parameters

WSL (

a

disp disp ind ∆U SL ∆U SL T ∆S SL

mJ ) m2

(

mJ ) m2

(

mJ ) m2

(

mJ ) m2

I

ε CO = 0.1072, σ CO = 3.3215, αC = 0

93.6

127.2

33.6

0

II

ε CO = 0.1072, σ CO = 3.3215, αC = 1.139, aC = 1.507

97.5

127.2

33.6

III

ε CO = 0.1072, σ CO = 3.3215, 104.4 αC = 1.80, aC = 0.50

127.2

IV

ε CO = 0.1448, σ CO = 3.2148, 127.2 αC = 0

170.0

ind T ∆S SL

(

mJ ) m2

T ∆S SL θ ∆U SL (°)

0

0.26

72.5

8.3

4.4

0.28

69.3

33.6

41.8

31.0

0.38

63.2

42.8

0

0

0.25

40.0

The units of ε CO , σ CO , and α C are kcal/mol, Å , and Å3 , respectively, while aC is

dimensionless. The statistical uncertainty in the calculated values of WSL is less than 1%. The reason behind the large difference in the values of the contact angle θ and WSL obtained in Cases III and IV can be attributed to the interfacial entropy of water. Indeed, for Case III, T ∆SSL , which is the interfacial entropic contribution to WSL , is 38% of ∆U SL , which is the energetic 30 ACS Paragon Plus Environment

Page 31 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


contribution to WSL . On the other hand, T ∆SSL is only 25% of ∆U SL in Case IV. Another interesting observation can be made by comparing the total energetic contribution, ∆U SL , to WSL ind in Cases III and IV. Recall that ∆U SL = ∆U Sdisp L + ∆U SL . Table 1 shows that ∆U SL is equal to

169 mJ/m 2 and 170 mJ/m 2 for Cases III and IV, respectively, which indicates that the energetic contributions to WSL are similar in both cases, and not sensitive to whether polarization effects are modeled explicitly (as in Case III) or implicitly (as in Case IV). Because, in Case III, the interfacial entropic contribution largely cancels the energetic contribution to WSL , the calculated value of WSL in Case III is lower than that in Case IV, and consequently, the contact angle is higher in Case III than in Case IV. Our result for the water contact angle in Case III (63.2 ° ), where polarization effects were modeled self-consistently, is in excellent agreement with the contact angle value reported by Kozbial et al. (65° ± 2°).19 In fact, even the values of the water contact angles obtained in Case I ( α C = 0 ) and Case II ( α C and aC were fitted to reproduce the polarizability tensor of graphene) (see Table 1) are in relatively good agreement with the recently reported experimental contact angle data.15–19 This follows because when polarization effects are modeled self-consistently (Cases II and III), the increase in the energetic component of the work of adhesion is cancelled, to a large extent, by the entropic component of the work of adhesion. However, modeling the induction energy implicitly using a pair-wise additive LJ potential (Case IV) clearly leads to a significant underprediction of the water contact angle. Therefore, we present a more complex picture of the water-graphene interface, where the contact angle of water on graphitic surfaces is determined not by the magnitude of the total binding energy of a water molecule interacting with



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 45

graphene, but instead, by the precise decomposition of this binding energy into dispersion and polarization components. Furthermore, comparing the results obtained for Cases I and IV (see Table 1), we find that T ∆S SL is not affected when polarization effects are not considered, although there is a large ∆U SL change in the calculated WSL values and the corresponding contact angles (the difference is ≈ 32° ) for these two cases. This again clearly shows that dispersion and polarization effects contribute differently to the wetting of graphitic surfaces by water. Physically, the interfacial entropic contribution to the work of adhesion represents the gain in entropy resulting from the reduced ordering of the interfacial water molecules when the interaction potential between the water molecules and graphite is gradually turned off. The large difference in the interfacial entropic component of WSL between Cases III and IV seems to indicate that polarization effects have a more pronounced effect on the structure of the interfacial water molecules compared to dispersion effects. Because the polarization of graphene results in a pronounced effect on the interfacial entropy of water, we next estimate the electric field exerted by the water molecules on the carbon atoms in graphene. Note that the magnitude of this electric field determines the magnitude of the induced dipole moments of the carbon atoms in graphene, which in turn, determines the magnitude of the interaction between the induced dipoles in graphene and the permanent multipoles in the water molecules. The contour plots for the x, y, and z components of the electric field exerted by the water molecules on the carbon atoms of the top-most graphene layer in graphite are shown in Figure 4. Note that these contour plots were calculated using MD simulations carried out using the LJ parameters, ε co = 0.1072 kcal/mol and σ CO = 3.3215 Å . The simulation methodology, including all the relevant simulation parameters, are identical to those 32 ACS Paragon Plus Environment

Page 33 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


used to calculate WSL using the FEP method. Note that the z direction is normal to the plane of the graphene layer. Because the graphene layer is periodic in the xy plane, the deduced values of the x and y components of the electric field are very similar, while the z component is considerably larger than both. Figure 4 shows that the x, y, and z components of the electric field progressively increase with an increase in the static dipole polarizability of the carbon atoms in graphene. Indeed, by carrying out further analysis, we find that the electric field depends exponentially on α C (see Supporting Information, Section S6 for additional details). Moreover, we also calculated the electric field felt by the carbon atoms in the second and third graphene layers, and the contour plots for the electric fields, calculated using the polarizability parameters,

α C = 1.80 Å3 and aC = 0.50 , are shown in Figure 5. Comparing Figure 4 and Figure 5, we find that the electric field felt by the carbon atoms in the second and third graphene layers is about

and

1 4

1 , respectively, of the electric field felt by the carbon atoms in the first graphene layer. 10

Because U pol ∝ E 2 ,20 our results show that the polarization of the carbon atoms in the first graphene layer makes the dominant contribution to the induction energy of the system.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 45

Figure 4. Contour plots showing the calculated electric field exerted by the water molecules on the carbon atoms in the top-most graphene layer. Images of the carbon atoms of the graphene layer are superimposed on the contour plots for better visualization. The rows correspond to different values of the polarizability parameters, αc and ac , while the columns correspond to the

x, y, and z components of the electric field vector. The units of αC and components of the electric field vector are Å3 and MV/m, respectively, while aC is dimensionless. The graphene layer measures 27 Å and 25.5 Å in the x and y directions, respectively. The mean and standard deviation of the electric field are indicated on the top of each of the plots.


Page 35 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Electric field exerted by the SWM4-NDP water molecules on the carbon atoms in graphite. In (a), we show the components of the electric field exerted on the second graphene layer, and in (b), we show the components of the electric field exerted on the third graphene layer. The averages as well as the first standard deviations of the components of the electric field over all the carbon atoms in the periodic graphene layer are indicated on top of each of the subplots in (a) and (b).

To precisely pinpoint the reason behind the dependence of the electric field on the static dipole polarizability of the carbon atoms in graphene, we calculated the induction energy of a single water molecule in various orientations interacting with graphene at a separation distance of 3.5 Å (see Figure 6 (a)). In Figure 6 (b), we also plotted the probability distribution of the different orientations of the water molecules in the first hydration layer (defined as the region until the first minimum in the density profile of the water molecules, located at around 5 Å from the top of the first graphene layer) next to graphite. An examination of Figure 6 (a) clearly shows that the minimum value of the induction energy corresponds to the 1 leg orientation of the water molecules, where one of the two OH bonds of water points towards the graphene surface. In turn, 35 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 45

this implies that the electric field felt by the carbon atoms in graphene is the highest for this particular orientation of the water molecule. Consequently, as shown in Figure 6 (b), upon progressively increasing the value of αC , the fraction of water molecules in the 1 leg orientation increases (this finding is in qualitative agreement with a previous study by Ho and Striolo22), because the water molecules at the interface seek to minimize their induction energy by adopting this particular orientation. However, this also results in an enhancement in the ordering of the interfacial water molecules, and hence, in a corresponding loss in the configurational degrees of freedom of the water molecules. Consequently, there is a pronounced loss of interfacial entropy upon increasing the value of αC . On the other hand, Figure 6 (b) shows that increasing the strength of the LJ interaction has no effect on the orientation of the interfacial water molecules, thereby explaining why

T ∆S SL is not affected when polarization effects are not considered (see ∆ U SL

Cases I and IV in Table 1). The results presented in this study are further corroborated by contact angle simulations of a cylindrical water nanodroplet on graphite (see Supporting Information, Sections S7 and S8 for additional details). The simulated contact angle corresponding to Case IV in Table 1 was found to be 18.5° lower than the one corresponding to Case III, thereby further validating our conclusion that polarization and dispersion effects contribute differently to the wetting of graphitic surfaces by water.


Page 37 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. (a) Calculated induction energy of a water molecule interacting with graphene as a function of cosψ , where ψ is the angle that the –OH bond of the water molecule makes with the normal to the plane of graphene. (b) Calculated probability of water molecules in different orientations in the first hydration layer next to graphite. The units of ε co , σ CO , and αC are

kcal/mol, Å , and Å3 , respectively. In the inset, we show a post-equilibrium MD simulation snapshot of 1000 SWM4-NDP water molecules interacting with graphite.

CONCLUSIONS

In summary, our study shows that when a film of water interacts with graphite, the polarization of graphene resulting from the electric field exerted by the water molecules has a pronounced effect on the interfacial entropy of water, an effect which is not observed when water molecules interact with graphene solely via dispersion interactions. Consequently, the contact angle of water on graphite deduced by self-consistently modeling the polarization of graphene using the 37 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 45

classical Drude oscillator model is remarkably different from that deduced by artificially augmenting the magnitude of the LJ potential to model the induction energy. Accordingly, our study highlights the necessity to self-consistently model polarization effects in order to make any quantitative predictions on the wetting properties of any polarizable material like graphene by a polar solvent like water. Moreover, by calculating the electric field exerted by the water molecules on graphene, we discovered a strong dependence of this electric field on the static dipole polarizability of the carbon atoms in graphene. We believe that such insights, which are challenging to reproduce experimentally, will also be useful in future fundamental studies on the role of polarization effects on the friction coefficient of water on graphene and inside carbon nanotubes.79 Furthermore, it would be interesting to extend the results and analyses presented in this study to other 2D materials like hexagonal boron nitride (hBN), transition metal dichalcogenides (e.g., MoS2), and phosphorene. For example, by carrying out ab initio simulations using DFT, Kumar et al.80 have recently reported that the polarizability of monolayer hexagonal boron nitride (hBN) is similar to that of graphene, while that of a phosphorus atom in phosphorene is three times larger than the polarizability of a carbon atom in graphene.81 We hope that our work will stimulate the use of polarizable force fields to describe the interaction of water molecules with other 2D materials like hBN, MoS2, and phosphorene, which will provide new insights for the use of these 2D materials in applications such as seawater desalination2 and osmotic power harvesting.82,83

ASSOCIATED CONTENT

Supporting Information

Further details on the computational setup used for carrying out the ab initio and force fieldbased MD simulations reported in this Article is included in the Supporting Information document. This material is available free of charge via the Internet at http://pubs.acs.org/. 38 ACS Paragon Plus Environment

Page 39 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


AUTHOR INFORMATION

Corresponding Author *E-mail: [email protected] NOTES The authors declare no competing financial interest

ACKNOWLEDGEMENTS

The authors thank Ananth G. Rajan, Zhe Yuan, Hari Katepalli, Vishnu Sresht, Annalisa Cardellini, Yamini Krishnan and Agilio Padua for insightful discussions. This work was supported by the National Science Foundation (NSF) grant number CBET-1511526. For computational resources, this work utilized the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF grant number ACI-1548562.

REFERENCES

(1)

Cohen-Tanugi, D.; Grossman, J. C. Water Desalination across Nanoporous Graphene. Nano Lett. 2012, 12 (7), 3602–3608.

(2)

Surwade, S. P.; Smirnov, S. N.; Vlassiouk, I. V.; Unocic, R. R.; Veith, G. M.; Dai, S.; Mahurin, S. M. Water Desalination Using Nanoporous Single-Layer Graphene. Nat. Nanotechnol. 2015, 10 (5), 459–464.

(3)

Dhiman, P.; Yavari, F.; Mi, X.; Gullapalli, H.; Shi, Y.; Ajayan, P. M.; Koratkar, N. Harvesting Energy from Water Flow over Graphene. Nano Lett. 2011, 11 (8), 3123–3127.

(4)

Garaj, S.; Hubbard, W.; Reina, A.; Kong, J.; Branton, D.; Golovchenko, J. A. Graphene as a Subnanometre Trans-Electrode Membrane. Nature 2010, 467 (7312), 190–193.

(5)

Schneider, G. F.; Kowalczyk, S. W.; Calado, V. E.; Pandraud, G.; Zandbergen, H. W.; Vandersypen, L. M. K.; Dekker, C. DNA Translocation through Graphene Nanopores. Nano Lett. 2010, 10 (8), 3163–3167.

(6)

Merchant, C. A.; Healy, K.; Wanunu, M.; Ray, V.; Peterman, N.; Bartel, J.; Fischbein, M. D.; Venta, K.; Luo, Z.; Johnson, A. T. C.; et al. DNA Translocation through Graphene Nanopores. Nano Lett. 2010, 10 (8), 2915–2921.

(7)

Algara-Siller, G.; Lehtinen, O.; Wang, F. C.; Nair, R. R.; Kaiser, U.; Wu, H. A; Geim, A K.; Grigorieva, I. V. Square Ice in Graphene Nanocapillaries. Nature 2015, 519 (7544), 443–445.

(8)

Chen, J.; Schusteritsch, G.; Pickard, C. J.; Salzmann, C. G.; Michaelides, A. Two Dimensional Ice from First Principles: Structures and Phase Transitions. Phys. Rev. Lett. 2016, 116 (2), 25501. 39 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 45

(9)

Rafiee, J.; Mi, X.; Gullapalli, H.; Thomas, A. V.; Yavari, F.; Shi, Y.; Ajayan, P. M.; Koratkar, N. A. Wetting Transparency of Graphene. Nat. Mater. 2012, 11 (3), 217–222.

(10)

Shih, C.-J.; Wang, Q. H.; Lin, S.; Park, K.-C.; Jin, Z.; Strano, M. S.; Blankschtein, D. Breakdown in the Wetting Transparency of Graphene. Phys. Rev. Lett. 2012, 109 (17), 176101.

(11)

Shih, C.-J.; Strano, M. S.; Blankschtein, D. Wetting Translucency of Graphene. Nat. Mater. 2013, 12 (10), 866–869.

(12)

Raj, R.; Maroo, S. C.; Wang, E. N. Wettability of Graphene. Nano Lett. 2013, 13 (4), 1509–1515.

(13)

Fowkes, F. M.; Harkins, W. D. The State of Monolayers Adsorbed at the Interface Solid— Aqueous Solution. J. Am. Chem. Soc. 1940, 62 (12), 3377–3386.

(14)

Morcos, I. On Contact Angle and Dispersion Energy of the Cleavage Graphite/water System. J. Colloid Interface Sci. 1970, 34 (3), 469–471.

(15)

Li, Z.; Wang, Y.; Kozbial, A.; Shenoy, G.; Zhou, F.; McGinley, R.; Ireland, P.; Morganstein, B.; Kunkel, A.; Surwade, S. P.; et al. Effect of Airborne Contaminants on the Wettability of Supported Graphene and Graphite. Nat. Mater. 2013, 12 (10), 925–931.

(16)

Amadei, C. A.; Lai, C.-Y.; Heskes, D.; Chiesa, M. Time Dependent Wettability of Graphite upon Ambient Exposure: The Role of Water Adsorption. J. Chem. Phys. 2014, 141 (8), 84709.

(17)

Wei, Y.; Jia, C. Q. Intrinsic Wettability of Graphitic Carbon. Carbon N. Y. 2015, 87 (C), 10–17.

(18)

Ondarçuhu, T.; Thomas, V.; Nuñez, M.; Dujardin, E.; Rahman, A.; Black, C. T.; Checco, A. Wettability of Partially Suspended Graphene. Sci. Rep. 2016, 6 (1), 24237.

(19)

Kozbial, A.; Trouba, C.; Liu, H.; Li, L. Characterization of the Intrinsic Water Wettability of Graphite Using Contact Angle Measurements: Effect of Defects on Static and Dynamic Contact Angles. Langmuir 2017, 33 (4), 959–967.

(20)

Stone, A. The Theory of Intermolecular Forces; Oxford University Press: Oxford, U.K., 2013.

(21)

Heinz, H.; Jha, K. C.; Luettmer-Strathmann, J.; Farmer, B. L.; Naik, R. R. Polarization at Metal-Biomolecular Interfaces in Solution. J. R. Soc. Interface 2011, 8 (55), 220–232.

(22)

Ho, T. A.; Striolo, A. Polarizability Effects in Molecular Dynamics Simulations of the Graphene-Water Interface. J. Chem. Phys. 2013, 138 (5), 54117.

(23)

Hughes, Z. E.; Tomásio, S. M.; Walsh, T. R. Efficient Simulations of the Aqueous BioInterface of Graphitic Nanostructures with a Polarisable Model. Nanoscale 2014, 6 (10), 5438–5448.

(24)

Tocci, G.; Joly, L.; Michaelides, A. Friction of Water on Graphene and Hexagonal Boron Nitride from Ab Initio Methods: Very Different Slippage Despite Very Similar Interface Structures. Nano Lett. 2014, 14 (12), 6872–6877. 40 ACS Paragon Plus Environment

Page 41 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(25)

Li, H.; Zeng, X. C. Wetting and Interfacial Properties of Water Nanodroplets in Contact with Graphene and Monolayer Boron–Nitride Sheets. ACS Nano 2012, 6 (3), 2401–2409.

(26)

Willow, S. Y.; Zeng, X. C.; Xantheas, S. S.; Kim, K. S.; Hirata, S. Why Is MP2-Water “Cooler” and “Denser” than DFT-Water? J. Phys. Chem. Lett. 2016, 7 (4), 680–684.

(27)

Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a LongRange Dispersion Correction. J. Comput. Chem. 2006, 27 (15), 1787–1799.

(28)

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 (15), 154104.

(29)

Nagata, Y.; Ohto, T.; Bonn, M.; Kühne, T. D. Surface Tension of Ab Initio Liquid Water at the Water-Air Interface. J. Chem. Phys. 2016, 144 (20), 204705.

(30)

Leroy, F.; Liu, S.; Zhang, J. Parametrizing Nonbonded Interactions from Wetting Experiments via the Work of Adhesion: Example of Water on Graphene Surfaces. J. Phys. Chem. C 2015, 119 (51), 28470–28481.

(31)

Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 (2), 926–935.

(32)

Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123 (23), 234505.

(33)

Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press:New York, 2002.

(34)

Leroy, F.; Müller-Plathe, F. Dry-Surface Simulation Method for the Determination of the Work of Adhesion of Solid–Liquid Interfaces. Langmuir 2015, 31 (30), 8335–8345.

(35)

Zwanzig, R. W. High‐Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954, 22 (8), 1420–1426.

(36)

Adam, N. K. Use of the Term “Young”s Equation’ for Contact Angles. Nature 1957, 180 (4590), 809–810.

(37)

Taherian, F.; Leroy, F.; van der Vegt, N. F. A. Interfacial Entropy of Water on Rigid Hydrophobic Surfaces. Langmuir 2013, 29 (31), 9807–9813.

(38)

Taherian, F.; Marcon, V.; van der Vegt, N. F. A.; Leroy, F. What Is the Contact Angle of Water on Graphene? Langmuir 2013, 29 (5), 1457–1465.

(39)

Lamoureux, G.; Harder, E.; Vorobyov, I. V.; Roux, B.; MacKerell, A. D. A Polarizable Model of Water for Molecular Dynamics Simulations of Biomolecules. Chem. Phys. Lett. 2006, 418 (1–3), 245–249.

(40)

Lamoureux, G.; Roux, B. Modeling Induced Polarization with Classical Drude Oscillators: Theory and Molecular Dynamics Simulation Algorithm. J. Chem. Phys. 2003, 119 (6), 3025–3039. 41 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 45

(41)

Thole, B. T. Molecular Polarizabilities Calculated with a Modified Dipole Interaction. Chem. Phys. 1981, 59 (3), 341–350.

(42)

Szalewicz, K. Symmetry-Adapted Perturbation Theory of Intermolecular Forces. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2 (2), 254–272.

(43)

Parker, T. M.; Burns, L. A.; Parrish, R. M.; Ryno, A. G.; Sherrill, C. D. Levels of Symmetry Adapted Perturbation Theory (SAPT). I. Efficiency and Performance for Interaction Energies. J. Chem. Phys. 2014, 140 (9), 94106.

(44)

Turney, J. M.; Simmonett, A. C.; Parrish, R. M.; Hohenstein, E. G.; Evangelista, F. A.; Fermann, J. T.; Mintz, B. J.; Burns, L. A.; Wilke, J. J.; Abrams, M. L.; et al. Psi4: An Open-Source Ab Initio Electronic Structure Program. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2 (4), 556–565.

(45)

Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117 (1), 1–19.

(46)

Li, S.; Li, Q.; Carpick, R. W.; Gumbsch, P.; Liu, X. Z.; Ding, X.; Sun, J.; Li, J. The Evolving Quality of Frictional Contact with Graphene. Nature 2016, 539 (7630), 541– 545.

(47)

Padua, A. A. H. fftool: A Tool to Build Force Field Input Files for Molecular Dynamics. https://github.com/ agiliopadua/fftool (accessed June 1, 2016).

(48) Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30 (13), 2157–2164. (49)

Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: A Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94 (26), 8897–8909.

(50)

Ma, M.; Tocci, G.; Michaelides, A.; Aeppli, G. Fast Diffusion of Water Nanodroplets on Graphene. Nat. Mater. 2015, 15 (1), 66–71.

(51)

Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31 (3), 1695–1697.

(52)

Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; CRC Press: Boca Raton, FL, 1988.

(53)

Wu, Y.; Aluru, N. R. Graphitic Carbon–Water Nonbonded Interaction Parameters. J. Phys. Chem. B 2013, 117 (29), 8802–8813.

(54)

Fitzner, M.; Joly, L.; Ma, M.; Sosso, G. C.; Zen, A.; Michaelides, A. Communication: Truncated Non-Bonded Potentials Can Yield Unphysical Behavior in Molecular Dynamics Simulations of Interfaces. J. Chem. Phys. 2017, 147 (12), 121102.

(55)

Jenness, G. R.; Jordan, K. D. DF-DFT-SAPT Investigation of the Interaction of a Water Molecule to Coronene and Dodecabenzocoronene: Implications for the Water - Graphite Interaction. J. Phys. Chem. C 2009, 113 (23), 10242–10248.

(56)

Jenness, G. R.; Karalti, O.; Jordan, K. D. Benchmark Calculations of Water–acene 42 ACS Paragon Plus Environment

Page 43 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Interaction Energies: Extrapolation to the Water–graphene Limit and Assessment of Dispersion–corrected DFT Methods. Phys. Chem. Chem. Phys. 2010, 12 (24), 6375. (57)

Ma, J.; Alfè, D.; Michaelides, A.; Wang, E. The Water-Benzene Interaction: Insight from Electronic Structure Theories. J. Chem. Phys. 2009, 130 (15), 154303.

(58)

Al-Hamdani, Y. S.; Alfè, D.; Michaelides, A. How Strongly Do Hydrogen and Water Molecules Stick to Carbon Nanomaterials? J. Chem. Phys. 2017, 146 (9), 94701.

(59)

Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77 (18), 3865–3868.

(60)

Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Van Der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92 (24), 246401.

(61)

Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Erratum: Van Der Waals Density Functional for General Geometries [Phys. Rev. Lett. 92 , 246401 (2004)]. Phys. Rev. Lett. 2005, 95 (10), 109902.

(62)

Lee, K.; Murray, É. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy van Der Waals Density Functional. Phys. Rev. B 2010, 82 (8), 81101.

(63)

Ma, J.; Michaelides, A.; Alfè, D.; Schimka, L.; Kresse, G.; Wang, E. Adsorption and Diffusion of Water on Graphene from First Principles. Phys. Rev. B 2011, 84 (3), 33402.

(64)

Stone, A. J. Distributed Multipole Analysis: Stability for Large Basis Sets. J. Chem. Theory Comput. 2005, 1 (6), 1128–1132.

(65)

Tafipolsky, M.; Engels, B. Accurate Intermolecular Potentials with Physically Grounded Electrostatics. J. Chem. Theory Comput. 2011, 7 (6), 1791–1803.

(66)

Singla, S.; Anim-Danso, E.; Islam, A. E.; Ngo, Y.; Kim, S. S.; Naik, R. R.; Dhinojwala, A. Insight on Structure of Water and Ice Next to Graphene Using Surface-Sensitive Spectroscopy. ACS Nano 2017, 11 (5), 4899–4906.

(67)

Yu, E. K.; Stewart, D. A.; Tiwari, S. Ab Initio Study of Polarizability and Induced Charge Densities in Multilayer Graphene Films. Phys. Rev. B 2008, 77 (19), 195406.

(68)

Phillips, J. Chemical Theory of Bonds and Energy Gaps In Covalent Bonding in Crystals, Molecules and Polymers; The University of Chicago Press:Chicago, 1969; pp 125-130

(69)

Benedict, X.; Louie, S. G.; Cohen, M. L. Static Polarizabilities of Single-Walled Carbon Nanotubes. Phys. Rev. B 1995, 52 (11), 8541–8549.

(70)

Kozinsky, B.; Marzari, N. Static Dielectric Properties of Carbon Nanotubes from First Principles. Phys. Rev. Lett. 2006, 96 (16), 166801.

(71)

Karapetian, K.; Jordan, K. D. Properties of Water Clusters on a Graphite Sheet. In Water in Confining Geometries; 2003; pp 139–150.

(72)

Langlet, R.; Devel, M.; Lambin, P. Computation of the Static Polarizabilities of MultiWall Carbon Nanotubes and Fullerites Using a Gaussian Regularized Point Dipole Interaction Model. Carbon N. Y. 2006, 44 (14), 2883–2895. 43 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 45

(73)

Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van Der Waals Complexes. Chem. Rev. 1994, 94 (7), 1887–1930.

(74)

Israelachvili, J. N. Intermolecular and Surface Forces; Academic Press Inc.:San Diego, CA, 2011

(75)

Dequidt, A.; Devémy, J.; Pádua, A. A. H. Thermalized Drude Oscillators with the LAMMPS Molecular Dynamics Simulator. J. Chem. Inf. Model. 2016, 56 (1), 260–268.

(76)

Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The Missing Term in Effective Pair Potentials. J. Phys. Chem. 1987, 91 (24), 6269–6271.

(77)

Yu, H.; Karplus, M. A Thermodynamic Analysis of Solvation. J. Chem. Phys. 1988, 89 (4), 2366–2379.

(78)

Vargaftik, N. B.; Volkov, B. N.; Voljak, L. D. International Tables of the Surface Tension of Water. J. Phys. Chem. Ref. Data 1983, 12 (3), 817–820.

(79)

Falk, K.; Sedlmeier, F.; Joly, L.; Netz, R. R.; Bocquet, L. Molecular Origin of Fast Water Transport in Carbon Nanotube Membranes: Superlubricity versus Curvature Dependent Friction. Nano Lett. 2010, 10 (10), 4067–4073.

(80)

Kumar, P.; Chauhan, Y. S.; Agarwal, A.; Bhowmick, S. Thickness and Stacking Dependent Polarizability and Dielectric Constant of Graphene–Hexagonal Boron Nitride Composite Stacks. J. Phys. Chem. C 2016, 120 (31), 17620–17626.

(81)

Kumar, P.; Bhadoria, B. S.; Kumar, S.; Bhowmick, S.; Chauhan, Y. S.; Agarwal, A. Thickness and Electric-Field-Dependent Polarizability and Dielectric Constant in Phosphorene. Phys. Rev. B 2016, 93 (19), 195428.

(82)

Siria, A.; Poncharal, P.; Biance, A.-L.; Fulcrand, R.; Blase, X.; Purcell, S. T.; Bocquet, L. Giant Osmotic Energy Conversion Measured in a Single Transmembrane Boron Nitride Nanotube. Nature 2013, 494 (7438), 455–458.

(83)

Feng, J.; Graf, M.; Liu, K.; Ovchinnikov, D.; Dumcenco, D.; Heiranian, M.; Nandigana, V.; Aluru, N. R.; Kis, A.; Radenovic, A. Single-Layer MoS2 Nanopores as Nanopower Generators. Nature 2016, 536 (7615), 197–200.


Page 45 of 45

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TOC Graphic


Insights on the Role of Many-Body Polarization Effects in the Wetting

Recommend Documents