Article Cite This: J. Phys. Chem. C 2017, 121, 28166−28179
pubs.acs.org/JPCC
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, United States
J. Phys. Chem. C 2017.121:28166-28179. Downloaded from pubs.acs.org by UNIV OF SOUTH DAKOTA on 09/03/18. For personal use only.
S Supporting Information *
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 out 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 that is remarkably different from the contact angle that 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 interaction. Consequently, polarization and dispersion interactions 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.
■
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 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, pairwise 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 constituting 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
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 and nanopore-based DNA sequencing4−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 that 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 84−86°. However, recent experimental data suggest that (1) graphite is mildly hydrophilic, with a water contact angle in the range 61−68° 15−18 and (2) the higher © 2017 American Chemical Society
Received: September 6, 2017 Revised: November 11, 2017 Published: November 12, 2017 28166
DOI: 10.1021/acs.jpcc.7b08891 J. Phys. Chem. C 2017, 121, 28166−28179
Article
The Journal of Physical Chemistry C
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, 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 fieldbased 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 fieldbased 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 that 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-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 selfconsistent treatment of polarization effects instead of using the
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 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). Though 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 air−water 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. 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 freeenergy 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) and γL is the surface tension of water, and θ is the contact angle of water on graphene. Note that the solid−liquid work of 28167
DOI: 10.1021/acs.jpcc.7b08891 J. Phys. Chem. C 2017, 121, 28166−28179
Article
The Journal of Physical Chemistry C
Accordingly, Parker et al.43 recommended using sSAPT0/jun-ccpVDZ as the bronze standard for SAPT 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 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 that allowed us to calculate the binding energy of a water molecule with larger PAHs like the dicircumcoronene molecule (C96H24). 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, 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 for additional details). 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 (1 K), 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
conventional approach in previous MD simulation studies where the interaction energy between water molecules and graphene has been modeled using a pairwise 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 the dispersion interaction, 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 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) E SAPT0 = [Eelst ]elst + [Eexch ]exch (20) (20) (2) + [E ind + Eexch − ind + δE HF ]ind (20) (20) + [Edisp + Eexch − disp]disp
(2)
SAPT0
where E is the binding energy that 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 Parker et al.43 Each of the terms in eq 2 is of the form Eij, 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 as sSAPT0, by moderately modifying the (20) SAPT0 E(20) exch−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-cc-pVDZ, Parker et al.43 reported obtaining very accurate binding energies using sSAPT0/jun-cc-pVDZ 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). 28168
DOI: 10.1021/acs.jpcc.7b08891 J. Phys. Chem. C 2017, 121, 28166−28179
Article
The Journal of Physical Chemistry C 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 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 identical to the one that we used to carry out the work of adhesion simulations. With 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.
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). 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 on the basis of the lateral position of the water molecule (difference in the total binding energy
