This is an open access article published under an ACS AuthorChoice License, which permits copying and redistribution of the article or any adaptations for non-commercial purposes.
Article pubs.acs.org/JPCB
Cite This: J. Phys. Chem. B 2018, 122, 3626−3634
Is Water at the Graphite Interface Vapor-like or Ice-like? Yuqing Qiu, Laura Lupi, and Valeria Molinero* Department of Chemistry, The University of Utah, 315 South 1400 East, Salt Lake City, Utah 84112-0850, United States
J. Phys. Chem. B 2018.122:3626-3634. Downloaded from pubs.acs.org by 191.101.54.50 on 09/03/18. For personal use only.
S Supporting Information *
ABSTRACT: Graphitic surfaces are the main component of soot, a major constituent of atmospheric aerosols. Experiments indicate that soots of different origins display a wide range of abilities to heterogeneously nucleate ice. The ability of pure graphite to nucleate ice in experiments, however, seems to be almost negligible. Nevertheless, molecular simulations with the monatomic water model mW with water-carbon interactions parameterized to reproduce the experimental contact angle of water on graphite predict that pure graphite nucleates ice. According to classical nucleation theory, the ability of a surface to nucleate ice is controlled by the binding free energy between ice immersed in liquid water and the surface. To establish whether the discrepancy in freezing efficiencies of graphite in mW simulations and experiments arises from the coarse resolution of the model or can be fixed by reparameterization, it is important to elucidate the contributions of the water− graphite, water−ice, and ice−water interfaces to the free energy, enthalpy, and entropy of binding for both water and the model. Here we use thermodynamic analysis and free energy calculations to determine these interfacial properties. We demonstrate that liquid water at the graphite interface is not ice-like or vapor-like: it has similar free energy, entropy, and enthalpy as water in the bulk. The thermodynamics of the water−graphite interface is well reproduced by the mW model. We find that the entropy of binding between graphite and ice is positive and dominated, in both experiments and simulations, by the favorable entropy of reducing the ice−water interface. Our analysis indicates that the discrepancy in freezing efficiencies of graphite in experiments and the simulations with mW arises from the inability of the model to simultaneously reproduce the contact angle of liquid water on graphite and the free energy of the ice−graphite interface. This transferability issue is intrinsic to the resolution of the model, and arises from its lack of rotational degrees of freedom. the other hand, cannot template the ice structure.5 Nevertheless, graphitic surfaces nucleate ice in simulations with the mW model.5−7,9,14−18 The simulations indicate that the ice nucleation rate as a function of the water−carbon attraction is nonmonotonous: it first increases, and then decreases sharply before further increasing again.15 This implies that the binding free energy of ice to graphite changes nonmonotonically as the graphitic surface becomes more hydrophilic. The second goal of this work is to compute the binding free energy of graphite to mW ice, determine its enthalpic and entropic contributions, and elucidate how these contributions change with the strength of the water−carbon attraction. The ability of graphite to nucleate ice in experiments, however, is not well established. Most evidence suggests that graphite has little or no ice nucleation ability,19−22 although soot containing graphitic surfaces can nucleate ice at temperatures up to 16 K above the homogeneous limit.19,23−25 These results suggest that the binding free energy of ice to graphite is closer to zero in experiments than in simulations. The third goal of this study is to elucidate which contribution to the free energy, entropy or enthalpy, is responsible for this difference between the experiments and the mW model, and to determine
1. INTRODUCTION It has been argued that hydrophobic surfaces promote vaporlike fluctuations in the interfacial water.1 Graphite is considered a prototypical example of hydrophobic surface,2,3 although some studies argue that it is more hydrophilic than previously expected.4 Molecular simulations with the monatomic water model mW indicate that the first layer of interfacial water on graphite has domains with order intermediate between ice and liquid,5 that these domains are the birthplace of ice nuclei,6 and that the degree of ice-like order at the surface increases with supercooling.7 This poses the question of whether water at the graphite interface is more or less ordered than in the bulk, and to which extent the simulations reproduce the experimental surface entropy of water at the graphite interface and its temperature dependence. Answering these questions is the first goal of this work. Classical nucleation theory (CNT) indicates that surfaces that nucleate ice must bind to ice.8 The stronger the ice nucleus binds to the surface, the lower is the free energy barrier for ice nucleation, and the warmer the temperature at which the surface promotes the nucleation of ice.8,9 The best ice nucleating agents, such as ice nucleating proteins10−12 and monolayers of long-chain alcohols,8,13 achieve a strong binding free energy to ice by providing ordered arrangements of hydroxyl groups that match the order of water molecules in specific faces of ice.8,10−13 The regular arrangement of OH groups results in strong hydrogen bonds between ice and the surface that acts as a template for ice formation. Graphite, on © 2018 American Chemical Society
Special Issue: Benjamin Widom Festschrift Received: November 21, 2017 Revised: January 2, 2018 Published: January 3, 2018 3626
DOI: 10.1021/acs.jpcb.7b11476 J. Phys. Chem. B 2018, 122, 3626−3634
Article
The Journal of Physical Chemistry B
plane, by mold integration using a Monte Carlo code.47 This validates our implementation of mold integration through molecular dynamics in LAMMPS. To compute the binding of ice to graphite at 273 K, we build a mold surface of the basal plane of ice composed of 288 potential energy wells with area 5 nm × 5 nm in the first contact layer on graphite. The mold and the graphite interface are exposed to bulk water containing 4116 molecules. The integral computed is γmold (r0w) = γi−w + γg−i − γg−w, where γi−w, γg−i, and γg−w are the surface tension of the ice−water, graphite−ice, and graphite−water interfaces, respectively. The maximum depth is set to εm = 8 RT; the averaged number of filled wells is computed with the same procedure as described for the ice−water interface. We find that r0w increases with increasing εw−c. We determine that for εw−c = 0.13 kcal mol−1, r0w = 0.48 ± 0.01 Å, and the free energy is integrated for rw = 1.00, 1.05, 1.10, 1.15, 1.20 Å. For εw−c = 0.20 kcal mol−1 we determine r0w = 0.77 ± 0.01 Å and integrate the free energy for rw = 1.05, 1.10, 1.15, 1.20, 1.25 Å. For εw−c = 0.25 kcal mol−1, we determine r0w = 0.88 ± 0.01 Å, and integrate the free energy for rw = 1.13, 1.15, 1.20, 1.22, 1.25 Å. Further increase of εw−c would result in a large r0w that would lead to the capture of more than one water molecule per well in the thermodynamic integration, invalidating the method. Hence, the mold integration cannot be used to compute the surface tension of strongly attractive graphitic surfaces. Figure S1 shows the fit and values of the γmold (r0w) for each εw−c. We compute the binding free energy ΔGbind of each surface from γmold (r0w) of that surface and the ice−liquid surface tension γi−w, ΔGbind = γg−i − γg−w − γi−w = γmold (r0w) − 2γi−w. We quantify the ice-like tetrahedral order of water and ice with the bond order parameter q3, as in ref 7. The average values of q3 for interfacial and bulk water as a function of temperatures are taken from ref 7. Here we compute the distribution of q3 of mW ice at the graphitic surfaces and compare it with that of bulk ice in order to estimate the surface entropy of the graphite−ice interface. We first crystallize a box of 4116 water molecules in contact with a 5 nm × 5 nm graphite surface. We select the water molecules in two bilayer of ice on each side of the graphitic surface, and compute their q3 values every 10 ps for 10 ns at 270 K. As the q3 order parameter considers the four closest neighbors to a given molecule, and the closest layer of ice to graphite is under coordinated, to make a fair comparison to the order parameter and its fluctuations at the graphite interface we select only two bilayers of ice from a trajectory of bulk ice, and compute their q3 value every 10 ps for 10 ns at 270 K, while only accounting for neighbors within the two bilayers. The resulting histograms of q3 for ice at the graphitic surfaces and in bulk are plotted in Figure S2. To measure the contact angle θ of mW water droplets on graphite and its dependence on temperature, we place a droplet containing 5241 water molecules on a 24 nm × 24 nm graphite surface. After equilibrating the droplet for 5 ns, we compute θ every 1 ps and averaged over 3 ns. The contact angles are computed as in ref 7, setting the temperature of the droplet to 218, 220, 230, 240, 250, 260, 270, 280, and 298 K.
whether the discrepancy arises from the interaction potentials or from the low resolution of the coarse-grained model.
2. SIMULATION MODELS AND METHODS Water is modeled with the monatomic water model mW,26 which has been extensively used for simulations of water and ice.5−8,14,26−44 Graphite is modeled as a single rigid layer of graphene following ref 5, and its interaction with water is represented by the two-body term of the Stillinger−Weber (SW) potential45 with water−carbon interaction size, σw−c = 0.32 nm, which reproduces the water−carbon distances of atomistic models.46 The strength of the water−carbon interaction, εw−c = 0.13 kcal mol−1, was parametrized5 to reproduce the contact angle of water droplets on graphite, θ = 86°.3 The other parameters in the SW potential for water− carbon are the same as in ref 45. We also consider graphitic surfaces with stronger water−carbon attractions, εw−c = 0.20 and 0.25 kcal mol−1, for which we compute the freezing temperatures following the same protocols as in ref 7. We compute the surface tensions of graphite−ice and ice− water interfaces with the mold integration method.47,48 To compute the surface tension of the ice−water interface, γice−water, we build a layer of 288 potential energy wells (“molds”) in a 5 nm × 5 nm area, where each mold is in the position of the oxygen in a single bilayer of the basal plane of ice. This mold is immersed in a simulation cell with dimensions 5 nm × 5 nm × 5 nm containing 4608 water molecules in the liquid state. We implement the interaction of the potential well and water as a continuous function, following refs 47,48, using the tabulated pair-potential function of LAMMPS.49 To determine the largest well radius r0w at which no induction time is required for the growth of ice, which we identify with the CHILL+ algorithm,34 we run simulations with various well radii, scanning from 0 to 1 Å, every 0.01 Å, with the depth of the well fixed at εm = 8 RT. For each radius we perform five independent molecular dynamics simulations starting with different random velocities. Each simulation is evolved for 1 ns in the NVT ensemble at 273 K, the temperature of ice−liquid coexistence for the mW model.35 The temperature is controlled with Nosé−Hoover thermostat50,51 with a relaxation time constant of 0.5 ps. The equations of motion are integrated with the velocity Verlet algorithm using periodic boundary conditions and a time step of 0.5 fs (this small step is required to integrate the steep interactions of mW with the mold). We find r0w = 0.53 ± 0.01 Å, smaller than the r0w = 1.00 Å reported for mold integration of the basal plane of mW using Monte Carlo simulations.47 We propagate the error bar in the surface tension from the error bar of r0w. The free energy of ordering water into ice through the molds is integrated for the potential wells width rw = 0.8, 0.9, 1.0, 1.1, 1.2 Å for ice−water interface, in each case turning on the potential energy wells from zero to its maximum depth, εm = 8 RT. The values of εm scanned for each system are listed in Table S1. For each value of εm of each rw, we count the number of wells filled with water every 1 ps and averaged this count over a 1 ns simulation. We verify that with the depth of εm = 8 RT, the averaged number of filled wells equals the number of wells in the system, 288. The free energy for each rw corresponds to twice the ice−water surface free energy, γi−w (rw). These values are then fitted to a straight line (Figure S1) to extrapolate to the thermodynamic value of γi−w = 35.0 ± 0.2 mJ m2 at r0w = 0.53 Å. This value agrees well with the γi−w = 34.5 ± 0.8 mJ m−2 determined for mW exposing also the basal
3. RESULTS AND DISCUSSION 3.1. Water at the Graphite Surface Has Similar Entropy and Enthalpy as in Bulk Water. The surface entropy of the graphite−water interface, Sg−w = (∂S /∂Ag−w)p,T = −(∂γg−w/∂T)p,Ag−w, is a measure of the difference in entropy 3627
DOI: 10.1021/acs.jpcb.7b11476 J. Phys. Chem. B 2018, 122, 3626−3634
Article
The Journal of Physical Chemistry B
γw−v for the coarse-grained model at 273 K are in good agreement with the experimental values (Table 1). As a result, the water−graphite surface tension γg−w for the mW model (−4.7 mJ m−2) is essentially the same as for water (−5.3 mJ m−2). The negative sign of γg−w indicates that water at the graphite interface at 273 K is slightly more stable than in the bulk. We obtain the surface entropy of liquid water at the graphite interface, Sg−w, by deriving eq 3 with respect to temperature:
of water at the surface and in the bulk of the liquid. If water at the graphite interface is ice-like, the surface entropy would be negative; if it is vapor-like, it would be positive. In what follows we use thermodynamic relations to derive an expression of γg−w in terms of the surface tension of the water−vapor interface γw−v and the contact angle θ of liquid water on graphite and compute the surface entropy of water at the graphite interface from the temperature dependence of these properties. Figure 1a sketches the process of adhesion of a slab of liquid water to the graphite surface. The free energy per area
⎛ ∂cosθ ⎞ ⎟ Sg − w = −Sw − v × cos θ + γw − v ⎜ ⎝ ∂T ⎠ p , A
where Sw−v = −
associated with that process is the adhesion free energy, ΔGadhesion, which we express in terms of the surface tensions of the graphite−water interface γg−w, the graphite−vapor interface γg−v, and the water-vapor interface γw−v: (1)
We write ΔGadhesion as a function of γw−v and the contact angle θ assuming, as in ref 52, that Young’s equation is valid and that the reversible work of spreading vapor on graphite, γg−v, is zero: ΔGadhesion = γg − w − γw − v = −γw − v × (1 + cos θ )
(2)
We then rearrange eq 2 to express γg−w as a function of γw−v and θ: γg − w = −γw − v × cos θ
∂T
is the surface entropy of the p, A w−v
water−vapor interface. The contact angle of mW water on graphite is independent of temperature in the range from 218 to 298 K (Figure S2). That reduces eq 4 to Sg−w = −Sw−v × cos θ. We are not aware of experimental measurements of the contact angle of water on graphite as a function of temperature. Here we assume θ is independent of temperature, at least in the immediacy of the melting temperature, also in experiments. The surface entropy of the liquid−vapor interface, Sw−v, is always positive for both water and mW. The magnitude of this quantity, however, is underestimated by the mW model because it does not have rotational degrees of freedom54 (Table 1 and Figure 2). We find that while the tetrahedral ordering on interfacial water increases on cooling; it does it in a way that mimics the ordering in bulk liquid water7 (Figure 3). This explains why the surface entropy of water at the graphite surface is quite insensitive to temperature (Figure 2). Comparing the surface entropies of water at the graphite and vapor interfaces (Table 1), it is clear that not only the signs are opposite but also that the magnitude of Sg−w is very small, because θ of water on graphite is close to 90°. We conclude that water at the graphite interface is neither ice-like nor vapor-like; its free energy, entropy, and enthalpy (Table 1) are very similar to that of water in the bulk liquid. The analysis of this section indicates that water at the graphite interface is only slightly more stable than in the bulk. The stabilization of water at the surface is driven by a favorable, but also small, surface enthalpy (Table 1). The coarse-grained model represents well the stability of water at the graphite interface, γg−w, because mW reproduces well the water-vapor surface tension and, by construction, it matches the experimental contact angle of water on graphite. The entropic and enthalpic contributions to the surface free energy of water at the graphite interface are also well reproduced by the mW model, but only because a small value of cos θ (∼0.07 for graphite) disguises the large differences in surface entropy of the water−vapor interface in mW and real water (Table 1). mW would not reproduce well the surface entropy and enthalpy of water at surfaces that are very hydrophilic or very hydrophobic if the parametrization is performed to reproduce the contact angle of liquid water with the surface. Alternative
Figure 1. Processes of adhesion of water to graphite and binding of ice to graphite discussed in this study. (a) Process of binding a slab of water to a graphitic surface; we call the free energy per area of this process the adhesion free energy, ΔGadhesion. (b) Process of binding a graphite surface initially immersed in liquid water to ice; we call the free energy per area for this process the binding free energy, ΔGbind.
ΔGadhesion = γg − w − γg − v − γw − v
∂γw − v
( )
(4)
g−w
(3)
The hydrophobicity of the surface determines the sign of the surface free energy between graphite and water, γg−w. Hydrophobic surfaces (θ > 90°; cos θ < 0) destabilize water and have γsurface−w > 0, while hydrophilic surfaces (θ < 90°; cos θ > 0) stabilize water and have γsurface−w < 0. Graphite is borderline between hydrophobic and hydrophilic, with contact angle θ = 86°.3 The angle θ and water−vapor surface tension
Table 1. Thermodynamics of the Water−Vapor and Water−Graphite Interfaces at 273 K γw−v (mJ m−2) water mW a
a
76 68b
γg−w (mJ m−2)
θ c
86° 86°
Sw−v (μJ m−2 K−1)
−5.3 −4.7
d
141 69e
Sg−w (μJ m−2 K−1)
Hg−w (mJ m−2 K−1)
−9.8 −4.8
−8.0 −6.0
Ref 53. bRef 54. cRef 55. dRef 53. eRef 54. 3628
DOI: 10.1021/acs.jpcb.7b11476 J. Phys. Chem. B 2018, 122, 3626−3634
Article
The Journal of Physical Chemistry B
parametrizations based on the thermodynamics of the interfaces may not reproduce the contact angle of very hydrophobic or hydrophilic surfaces.56 Fully atomistic water models generally reproduce the surface tension of the water− vapor interface within 20% of the experimental value, and reproduce well the surface entropy of that interface.54 From that, we expect that fully atomistic models parametrized to reproduce the contact angle of water on graphite (or any other rigid surface), will properly capture the experimental thermodynamics of interfacial water. 3.2. Binding of Graphite to Ice Is Favored by Entropy. The entropy associated with the binding of graphite from water to ice (Figure 1b) is ΔS bind = Sg − i − Sg − w − Si − w
(5)
where Sg−w is the graphite−water interface derived in Section 3.1 and Sg−i and Si−w are the surface entropies of the graphite− water and ice−water interfaces, respectively. Replacing Sg−w with eq 4 we obtain ΔS bind = Sg − i + Sw − v cos θ − Si − w
(6)
As water molecules at the ice−graphite interface have order and fluctuations similar to that in bulk ice (Figure S3), we approximate Sg−i ≈ 0. Sg−w is also small (Section 3.1), therefore the surface entropy of the ice−water interface, Si−w, controls the sign and magnitude of the entropy of binding graphite to ice, ΔSbind. We obtain the surface entropy of the water−ice interface, Si−w, from the temperature dependence of the ice−water surface tension γi−w, using Turnbull’s heuristic relation,57 γi−w = λΔHmvm−2/3, where λ is a constant equal to 0.32 for water, ΔHm is the melting enthalpy of ice, and vm its molar volume. This relation is exact for hard spheres58 and a very good approximation for water.59 Using Turnbull’s relation, the expression of the ice−liquid surface entropy becomes
Figure 2. Temperature dependence of surface entropies and binding entropies of (a) water and (b) the mW model. Entropy of binding of graphite to ice (black), entropy of the water−vapor interface (blue), entropy of the graphite−water interface (orange), and entropy of the ice−water interface (red). Si−w decreases on supercooling, following the anomalous behavior of the heat capacity of supercooled liquid water. The lowering of the ice−liquid surface entropy is reflected almost symmetrically in the entropy of binding ice to graphite, ΔSbind.
⎛ ∂γ ⎞ Si − w = − ⎜ i − w ⎟ ⎝ ∂T ⎠ p ⎛ ∂ΔHm ⎞ −2/3 2 −5/3⎛ ∂vm ⎞ ⎟ ΔH ⎟v = − λ⎜ + λvm ⎜ ⎝ ∂T ⎠ p m ⎝ ∂T ⎠ p m 3 = −λvm−2/3ΔCp +
2 −2/3 λvm αΔHm 3
(7)
where ΔCp is the difference in heat capacity between ice and liquid and α is the thermal expansion coefficient of ice. The small value of α makes the second term negligible (this would still be the case if vm referred to the volume of liquid water). Replacing λ vm−2/3 = γi−w(Tm)/ΔHm(Tm), the surface entropy of
Figure 3. Temperature dependence of the tetrahedral ice-like ordering (measured by the bond-order parameter q3) for mW water in the first layer in contact with graphite (black circles) is the same as for the bulk liquid (red circles). Adapted from ref 7 with permission.
Table 2. Thermodynamics of the Ice−Water Interface, Ice Melting, and of Binding Graphite to Ice at the Corresponding Melting Temperatures: 273.15 K for Water,61 273 K for mW,35 252 K for TIP4P/2005,62 and 232 K for TIP4P62 ΔHm (kJ mol−1) water mW TIP4P TIP4P/2005
a
6.0 5.3b 4.39c 4.85c
γi−w (mJ m−2) d
32 35 27.2e 28.9e
ΔCp (J K−1 mol−1) a
40.6 7.2b 22.2c 29.7c
Si−w (μJ m−2 K−1)
ΔSbind (μJ m−2 K−1)
−TΔSbind (kJ mol−1 nm−2)
Hi−w (mJ m−2 K−1)
−215 −44 −137 −177
225 49 145f 185f
−37.0 −8.0 −20.3g −28.1g
−26.7 23.0 −4.6 −15.7
a Ref 62. bRef 44. cRef 62. dAverage value of refs 63,64. eRef 47. fSi−w = 112 μJ m−2 K−1 for TIP4P and Si−w = 113 μJ m−2 K−1 for TIP4P/2005 are estimated from γw−v(T) provided in ref 65. We assume the contact angle for water on a graphite surface is parametrized to reproduce the experimental contact angle, θ = 86° and computed their binding entropy ΔSbind using eq 6.
3629
DOI: 10.1021/acs.jpcb.7b11476 J. Phys. Chem. B 2018, 122, 3626−3634
Article
The Journal of Physical Chemistry B Table 3. Thermodynamics of Binding Graphite and Other Graphitic Surfaces to Ice at 273 K system
εw−c (kcal mol−1)
water mW
graphite 0.13 0.20 0.25
θa
ΔTf (K)
ΔGbind (kJ mol−1 nm−2)
−TΔSbind (kJ mol−1 nm−2)
ΔHbind (kJ mol−1 nm−2)
γg−i (mJ m−2)
86° 86° 0 0
∼0 12 ± 2b 19 ± 2 21 ± 2
>−4 −19.5 −27.0 −29.4
−37.0 −8.0 −8.4 >−10.8
>20.1 −2.1
