Cassie–Baxter and Wenzel States on a ... - ACS Publications

Jun 18, 2012 - The measured free energy barrier separating the Cassie−Baxter from the ... Finally, we found that the Cassie−Baxter/Wenzel transiti...
0 downloads 0 Views 3MB Size
Download PDF
Article pubs.acs.org/Langmuir

Cassie−Baxter and Wenzel States on a Nanostructured Surface: Phase Diagram, Metastabilities, and Transition Mechanism by Atomistic Free Energy Calculations Alberto Giacomello,† Simone Meloni,*,‡ Mauro Chinappi,§ and Carlo Massimo Casciola† †

Dipartimento di Ingegneria Meccanica e Aerospaziale, Università di Roma La Sapienza, Via Eudossiana 18, 00184 Roma, Italy School of Physics, Room 302 UCD-EMSC, University College Dublin, Belfield, Dublin 4, Ireland § Dipartimento di Fisica, Università di Roma La Sapienza, P. le A. Moro 5, 00185 Roma, Italy ‡

S Supporting Information *

ABSTRACT: In this work, we study the wetting of a surface decorated with one nanogroove by a bulk Lennard-Jones liquid at various temperatures and densities. We used atomistic simulations aimed at computing the free energy of the stable and metastable states of the system, as well as the intermediate states separating them. We found that the usual description in terms of Cassie−Baxter and Wenzel states is insufficient, as the system presents two states of the Cassie−Baxter type. These states are characterized by different curvatures of the meniscus. The measured free energy barrier separating the Cassie−Baxter from the Wenzel state (and vice versa) largely exceeds the thermal energy, attesting the existence of Cassie−Baxter/Wenzel metastabilities. Finally, we found that the Cassie−Baxter/Wenzel transition follows an asymmetric path, with the formation of a liquid finger on one side of the groove and a vapor bubble on the opposite side.

■

INTRODUCTION Hydrophobic materials can be turned into superhydrophobic ones if their surface is decorated with micro- or nanocorrugations. This phenomenon, known as “lotus effect” from the properties and structure of lotus leaves, results in liquid droplets with high contact angles (that can reach almost complete dewetting θ ∼180°1−3) with negligible associated hysteresis.4−7 Superhydrophobicity is due to the formation of a Cassie− Baxter (CB) state,8 i.e., a state in which the liquid does not penetrate into the hollows of the corrugated surface and, consequently, faces a composite interface consisting of both solid and vapor (see Figure 1A).9 Another state exists, in which the liquid is in contact with the entire exposed surface of the solid, that is called Wenzel (W) state11 (Figure 1B). For hydrophobic substrates, both states result in an increase of the contact angle, but only the CB state brings low contact angle hysteresis. Another interesting property shown by bulk liquids in the CB state under flow conditions is that they present a nonzero average velocity at the liquid/wall interface, which takes the name of slip velocity. This phenomenon is commonly characterized in terms of slip length Ls.12 It is nowadays clear that for smooth surfaces Ls can hardly exceed the nanometer range,13,14 while for patterned surfaces in CB state Ls on the micrometer scale have been observed.15 © 2012 American Chemical Society

Figure 1. (A) Cartoon of the Cassie−Baxter state in a drop (top) and a bulk liquid (bottom): the liquid does not wet the entire surface but remains suspended on top of the surface asperities. (B) Same as for the Wenzel state: the liquid is in contact with the asperities of the surface. θCB and θW denote the contact angles on a rough substrate when the system is in the Cassie−Baxter and Wenzel state, respectively.

The above characteristics make the CB state appealing for technological applications, and superhydrophobic materials attracted significant attention for their possible application in the production of nonwettable/self-cleaning textiles,16,17 low friction microdevices,18 and many others (see Zhang et al.19 for Received: May 5, 2012 Revised: June 15, 2012 Published: June 18, 2012 10764

dx.doi.org/10.1021/la3018453 | Langmuir 2012, 28, 10764−10772

Langmuir

Article

an extensive review of techniques for producing superhydrophobic surfaces and their possible technological applications). This has promoted applied research aimed at developing decorated surfaces having a high contact angle, with an associated low hysteresis, a large slip length, and that can maintain the system in the CB state for the widest interval of external conditions. Fundamental research is more focused on the identification of the factors controlling the stability of the CB vs W states,20−23 a subject that has been investigated also by continuum24 and atomistic25,26 simulations. Another objective is to assess whether CB/W metastabilities can exist,17,27−29 i.e., whether a liquid can be kept in the CB state in conditions in which the thermodynamic stable state is W (and vice versa). This phenomenon can happen when the free energy surface presents one absolute minimum and one or more local minima separated from the former by large free energy barriers (i.e., as compared to the thermal energy kBT). Metastability can also have a technological importance, as in practice, it represents a way of extending the range of stability of the CB state. On the other hand, a negative consequence of metastability is that it might prevent/slow down the “dewetting” process (transition from W to CB). Thus, the identification of factors controlling the height of these barriers is a subject of intense research. Metastabilities and hysteresis make the investigation of the CB/W phase transition by experiments or (atomistic) simulations very difficult. A first step toward the complete understanding of this process would be the determination of the density (ρ)−temperature (T) equilibrium phase diagram.30 We aim at achieving this objective by computing the free energy of the CB vs W phases at various ρ and T by atomistic simulations. To this end, we will use the combined restrained molecular dynamics - parallel tempering31 (RMD-PT) approach, which is well-suited to dealing with systems bearing known and unknown metastabilities. In addition to the phase diagram, by computing the free energy along the path connecting the CB and W states, we quantify the relevant free energy barriers present in the system, clarifying whether CB/W metastabilities can exist in the investigated system. Finally, by analyzing the atomistic configurations along this path we try to identify the mechanism of the transition. Previous experiments and simulations have shown that the CB/W phase diagram and transition rate depend on many factors: the chemical nature of the liquid, the intrinsic hydrophobicity of the material composing the surface,32 whether the liquid is deposited as a nano/micro droplet or in bulk (see Figure 1), the characteristics of the corrugations2,33−36 of the surface (shape,37 width, and height25,26 of the corrugations), and many more. In the present work, we focus on fundamental aspects of the CB/W transition, in particular, on the ρ−T phase diagram of a groove of aspect ratio close to unit. We postpone to future studies the investigation of different aspect ratios and geometries on the thermodynamics and kinetics of the transition. We consider the case of a bulk Lennard-Jones liquid (LJ) contained between two planar walls with a nanoscopic decoration on one of them (see Figure 2). The surfaces, also formed by LJ particles, interact with the liquid via a modified LJ potential (see next section) that can be tuned to achieve the desired hydrophobicity. The essential ingredients giving rise to the hydrophobic effect are the liquid being close to phase coexistence with its vapor and liquid−solid interactions being significantly less attractive than liquid−liquid interactions.38 In

Figure 2. Cartoon of the simulation box. Data are reported in Lennard-Jones units. H denotes the distance between the walls, L the computational box width, and h and l the height and the width of the groove, respectively. The thickness of the simulation box is 11 and the groove goes all through it along the z direction. The dashed rectangle, of height cy and width cx, represents the cell used for the calculation of the collective variable (see text). The simulation box is periodically replicated along directions x, y, and z.

view of this, LJ liquids have often been employed in studying the hydrophobic effect. For this reason, and for its simplicity, we consider that it is suitable for a fundamental study. Anticipating our results, we mention that the CB/W picture is insufficient to represent our system as we identify two CB states, characterized by different values of the curvature of the meniscus. At intermediate T, the low curvature CB state is more stable, while the high curvature CB state becomes more stable at lower T. We show that, at a variance with what was predicted by macroscopic theories,39 at the nanoscale the CB state can also be stable in the case of moderate hydrophilic surfaces (θY ≲ 90°, with θY being Young contact angle). We provide evidence that metastable states separated by stable ones by large free energy barriers (10 kBT or more) are supported by the considered system over a wide range of temperatures. This confirms that rare event techniques are necessary to study the CB/W phase diagram/transition, as brute force simulations may remain trapped in these metastable states. Concerning the transition mechanism, we found that, at a variance with what was proposed in literature,21,40,41 the CB/W transition follows an asymmetric path, with the formation of a liquid finger on one side of the groove and a vapor bubble on the opposite side. This path can still be explained in terms of macroscopic theories once unnecessary assumptions are lifted.

■

METHODS

In this section, we give an overview of the restrained molecular dynamics - parallel tempering method (RMD-PT) method. The readers already acquainted with these tools, or those uninterested in the simulative details, may proceed directly to the Results and Discussion section. The free energy calculations are performed using the RMD-PT method (see Orlandini et al.31 for a detailed description). In RMD-PT, we combine a specialized version of the Maragliano and VandenEijnden TAMD42 method with parallel tempering (PT).43,44 Atoms are evolved according to the following equation of motion

mr ̈ = −∇r U k(r, z) + thermo(β)

(1) −1

where m is the physical mass, β = (kBT) , and thermo(β) indicates that the atoms are coupled to a thermostat at β. The potential Uk(r, z) = V(r) + k/2[ε(r) − z]2 is the sum of the physical (V(r)) and restraining potentials. ε(r) is a suitable collective variable (CV) and z is its target value. Let us consider the following average: 10765

dx.doi.org/10.1021/la3018453 | Langmuir 2012, 28, 10764−10772

Langmuir

Article

Figure 3. Main panel: F(z) profiles at selected temperatures and three densities: 0.742 (top), 0.732 (middle), and 0.725 (bottom). The free energy is reported in LJ units. The integration constant of the thermodynamic integration is set such that the free energy at z = 200 is equal to zero. In the three panels on the right, we report the enlargement of selected F(z) curves. The blue curve in the lower right diagram is shifted. For the sake of clarity, the error bars49 are reported only for the curves shown in the right-hand side panels. In general, the maximum error is ≤16.5, ≤ 12.5, and ≤15 for the curves at ρ = 0.742, 0.732, and 0.725, respectively.

=

swapping between the configurations of two replicas, denoted by χ and ν, is attempted and accepted with a probability

τ 1 dt k[ε(r(t )) − z] τ →∞ τ 0 − ∫ dr k(ε(r) − z) exp[− βU k(r, z)]

f k (z) = lim −

∫

Ak(z) ⎛ Ak(z) ⎞ ⎟ = −∇z β −1 ln⎜ ⎝ A ⎠

a(r χ → rν, rν → r χ ) = min{1, exp[− (βν − βχ )(U k(r χ , z) − U k(rν, z))]}

(3)

where βν = 1/kBTν and βχ = 1/kBTχ and rν and rχ are the configurations of the ν and χ replicas, respectively. Assuming that Tν ≫ Tχ, the replica ν has a higher probability to overcome (free) energy barriers possibly present in directions other than ε(r). This results in a higher efficiency of PT at sampling multimodal restrained probability density functions, thus solving the problem of unknown CVs bearing metastabilities. The RMD-PT simulations described above have been performed using the LAMMPS software package46 “driven” by the PLUMED code47 for rare event simulations. In LAMMPS, we implemented the modified LJ potential governing the liquid/walls interaction, and in PLUMED, we added a CV suitable to study the CB/W transition (see Results and Discussion). To make the calculation computationally more efficient, we used the linked-cell based reordering algorithm48 that, in the case of a liquid sample, allows us to reduce the computational time by up a factor of 4.

(2)

where A k(z) = ∫ dr exp[−βUk(r, z)] and A = ∫ dr exp[−βV(r)] is the canonical partition function of the real system. Since A is zindependent, its introduction does not affect our argument, but it is necessary for the interpretation of f k(z) as the derivative of the free energy. The second equality in eq 2 stems from the assumption that, apart from the direction ε(r), the system is ergodic. f k(z) is the derivative of Fk(z) = −β−1 ln[A k(z)/A ]. Noting that limk→∞ exp[−βk/2(ε(r)−z)2]/(2π/βk)1/2 = δ(ε(x) − z), in this limit we have f k(z) = ∇zFk(z) → −β−1∇z ln Pε(z). Here, Pε(z) = ∫ dr w(r)δ(ε(r) − z) is the probability that ε(r) = z, with w(r) denoting the canonical probability distribution. Recalling that the free energy of a CV is defined as F(z) = −β−1 ln Pε(z), we find that in the proper limits eq 2 is an estimate of the derivative of the free energy: ∇zF(z). The time integral in the second term of eq 2 can be estimated by MD, with atoms evolving according to eq 1. This proves that by restrained MD we can estimate the derivative of the free energy with respect to one CV. The free energy profile can then be computed by numerical integration of f k(z) (i.e., thermodynamic integration45). In the above reasoning, we made the critical assumption that the only CV bearing metastabilities is ε(r). When this is not the case, MD is not efficient at sampling the distribution exp[−βUk(r, z)]/A k(z). This problem can be solved by combining RMD with parallel tempering (PT).43,44 In PT, the simulated system is composed of M replicas of the original system. All replicas are evolved by RMD at the same z value but at different temperatures. From time to time, a

■

RESULTS AND DISCUSSION The system for which we study the CB/W phase diagram and transition consists of a LJ liquid contained between two planar walls shaped as shown in Figure 2. The walls are also made of LJ particles, which are kept fixed at the lattice positions of a face-centered cubic (FCC) crystal. The density of the walls’ LJ particles is higher than that of the corresponding crystal at the thermodynamic conditions of the system so as to prevent crystallization of the liquid. The liquid and walls particles 10766

dx.doi.org/10.1021/la3018453 | Langmuir 2012, 28, 10764−10772

Langmuir

Article

interact via a modified LJ potential of the type v(r) = 4εLJ[(σLJ/ r)12 − c(σLJ/r)6], where c is used to tune the “hydrophobicity” of the walls, εLJ defines a scale for the energy, σLJ for the length, and r is the interatomic distance. We set c = 0.7, corresponding to a Young contact angle θY∼85° (see Supporting Information). The simulation box is a parallelepiped of size 118 × 72 × 11 in LJ units (LJ units are used throughout the text) and the distance between the walls, i.e., the space available to the liquid, is 53. The groove where the CB/W transition takes place is ∼11 × 11 × 11, i.e., approximately 10 times the particle size, consistent with a nanocorrugation. To study the CB/W phase diagram/transition, we compute the free energy of the system as a function of a CV able to monitoring the phase of the system. In the following, we will denote by ε(r) this CV, where r is the 3N vector of the atomic positions. It is worth recalling that the free energy of a collective variable is nothing else but the probability to observe the corresponding state expressed in the logarithmic scale and in kBT units F(z) = −kBT ln Pε(z)

(4)

Figure 4. Probability density Pε(z) vs z at various temperatures and at the same densities of Figure 3. The color code used in this figure is the same as in Figure 3. The three phases identified in our simulations, that we named Cassie−Baxter 1, Cassie−Baxter 2, and Wenzel are highlighted in the bottom panel. The red arrow in the central panel points to a bimodal Pε(z), indicating that the corresponding temperature is the CB1/CB2 transition temperature at this density.

where F(z) is the free energy and kB is the Boltzmann constant. Here, Pε(z) = ∫ dr w(r)δ(ε(r) −z) is the probability that ε(r) = z, with w(r) denoting the invariant distribution of the atomistic system. In the present work, we consider a system in the canonical ensemble, and therefore, w(r) is the Boltzmann distribution. By combining RMD-PT with Kirkwood’s thermodynamic integration,45 we are able to compute the free energy as a function of z (see Methods). From the free energy profile, via eq 4, one can identify what is the most probable state in the given thermodynamic conditions. The CV used in this investigation counts the number of particles contained within a cell enclosing the groove (the dashed rectangle in Figure 2). This definition captures the main feature characterizing the CB and W states: the degree of filling of the groove. When the system is in the CB case, i.e., when the groove is empty, the value of this CV is low. On the contrary, when the system is in the W state the groove is filled by the liquid and the value of this CV is higher. By comparing the free energy at the z-values corresponding to these two states, we can establish whether, in the given external conditions (T and nominal density ρ50), the CB or W state is more probable. To implement this approach, we need to know the values of z corresponding to the relevant states of the system. These values, that depend on external conditions are, a priori, not known. This problem can be solved by computing the free energy over a wide z interval at each ρ−T point considered. Apparently, this approach requires computation of more data than necessary. However, we can take advantage of the knowledge of the free energy over a wide z-range to check whether these states correspond to (global/local) minima of the free energy and to compute the magnitude of the barrier separating them. This will enable us to establish whether the system could get trapped in a metastable state. Free Energy Profiles. In Figure 3, the free energy profiles are reported for various temperatures (T = 0.76−0.9) and densities (ρ = 0.725, 0.732, and 0.742). These ρ−T values belong to the liquid domain of a bulk LJ system.51 The corresponding probability density functions Pε(z), calculated from F(z) through eq 4, are shown in Figure 4. Error bars grow with z (see the right-hand panels of Figure 3), which is typical of quantities that are computed via integration. Nevertheless, a clear trend of the free energy curves emerges from our data. As

a general remark, we notice that the Wenzel state, corresponding to the minima of F(z) at high z values (see Figure 5), is favored by high densities and temperatures. This

Figure 5. Density field ρz(x) at z = 400 (left), z = 1000 (center), and z = 2000 (right). These density fields correspond to the two minima of the free energy in the top panel on the right-hand side of Figure 3 and to the minimum at z = 2000 of the central panel of the same figure.

result is consistent with previous experimental and computational findings.17,27,28 In particular, we notice that, at the highest density, it is impossible to obtain the CB state even at low T. In fact, at this ρ there are no minima of the free energy in the z interval corresponding to the CB state (see top panel of Figure 3). This seems to indicate that there exists a critical density above which the system is always in the “wet” state. At lower densities, the system can exist in three states. We identify the states at z = 200−500 and z = 800−1500 as two different CB states (denoted CB1 and CB2 in the following) and the state at z > 1900 as a W state. CB1 and CB2 are, indeed, two different thermodynamic states, as confirmed by the contemporary presence of two minima of the free energy in the CB range for some ρ−T values (see right-hand panels of Figure 3). The phase transition between these two states is made evident by looking at the Pε(z) vs ρ and T graph (Figure 4). The contemporary presence of multiple minima in F(z) at several ρ−T points is not reflected in a multimodal Pε(z) as the difference in the value of the free energy at these minima is much larger than the thermal energy, with the noticeable 10767

dx.doi.org/10.1021/la3018453 | Langmuir 2012, 28, 10764−10772

Langmuir

Article

contrary, in the CB2 state α decreases with T. The data also indicate that α increases with ρ. However, in this respect, they are less univocal, and we do not draw any definitive conclusion on this. Finally, from Figure 7 we deduce that α must approach 0° before the system undergoes the CB2/W transition. This finding is consistent with the Gibbs’ criterion for the advancement of the triple line beyond a sharp corner,53,54 which prescribes that the triple line remains pinned (with our definition of the tangent angle) for (ϕ − θY) ≤ α ≤ (180° − θY), where ϕ is the angle of the corner. In our case, with θY ≈ 90° and ϕ = 90°, this corresponds exactly to 0° ≤ α ≤ 90°. We now compare the reported wetting states with the available macroscopic models for contact angles of a liquid drop on a rough surface. If the drop is sufficiently large as compared to the characteristic scale of the roughness, it is possible to break down the problem to that of finding the macroscopic apparent contact angle and a “cell” problem solving the microscopic details of surface wetting.21,39 The latter task is one of finding the position of the liquid/vapor interface minimizing surface energy over a (periodic) cell on the scale of the roughness. It is possible to compare this microscopic problem with our system, that can be viewed as a model cell for a surface having nanoscopic (periodic) roughness. A rigorous result found by Alberti and De Simone39 is that close to θY ∼90°, on both the hydrophilic and hydrophobic sides, the thermodynamically stable state is the Wenzel one (in the cited derivation, metastable states are not explicitly discussed). At the nanoscale, this seems to be contradicted by our results, which show that, even for a slightly hydrophilic surface θY ≲ 90°, at least one state of the Cassie−Baxter type is thermodynamically stable over a wide range of temperatures (see Figure 4). Free Energy Barriers. Let us now discuss whether our results support the hypothesis, previously proposed in literature,2,27−29 of the existence of CB/W metastability. Metastabilities can occur when a large free energy barrier (ΔF†) separates the CB and W minima, at ρ−T points where they are both present. By large, we mean that ΔF† ≫ kBT. The relation between this barrier, computed from the data reported in Figure 3, and the one the system must overcome in the CB/ W transition deserves some clarification. In the present study, we use a CV that is able to identify the different states of the system but might be inadequate to describe the mechanism of the phase transition,55 a task that might require more CVs. This results in a free energy barrier that is lower than the real one (see Bolhuis et al.55 for a thorough discussion of this problem). Therefore, a large ΔF† estimated from our simulations is a very strong indication of the existence of the CB/W metastability. Figure 8 shows the free energy barrier ΔF† relative to the CB2 → W and the W → CB2 transitions as a function of ρ and T. The branches relevant to metastability are W → CB2 at low T, where the stable state is CB2, and CB2 → W at high T. Within these branches, ΔF† ranges between 10 and 30 kBT, with W → CB2 barriers typically higher than the CB2 → W ones. Such large ΔF† can explain the observed CB and W metastabilities. Moreover, the fact that the W → CB2 barriers are typically higher than the CB2 → W ones, of a factor 2 in some cases, might explain the experimental observation that the wetting process is, in several cases, irreversible.28,56−58 Mechanism of the CB/W Transition. It would be very useful to obtain information on the atomistic mechanism of the CB/W transition. For example, this might help develop strategies to widen the range of metastability of the CB2 state. Unfortunately, the mechanism is not readily available in

exception of ρ = 0.732/T = 0.77, highlighted in Figure 4 with a red arrow. The bimodality of Pε(z) indicates that, at this ρ and T, the CB1 and CB2 phases are in equilibrium. The presence of multiple minima of F(z) in the range z ∈ [800, 2000] supports the experimental and computational observations of CB/W metastability. Before discussing this point, in the next section we focus on the characteristics of the two CB and the W states. Characterization of the Wetting States. In Figure 5, we show the density field ρz(x)52 at z = 400, 1000, and 2000 corresponding to three minima in the top and central charts of the right-hand panel of Figure 3. It is apparent that the state at high z (∼2000) corresponds to a W state. The other two states are of the CB type. What distinguishes these states is the curvature of the meniscus. To make this argument more quantitative, we focus on the value of the angle formed by the tangent of the meniscus at the triple line with the wall, denoted by α in the following. α is ∼90° for CB1 and ranges between 0° and 50° for CB2. The dependency of α on ρ and T in an equilibrium sample is also of interest. This can be obtained from ρz(x) at z values corresponding to the maximum of Pε(z) at given ρ and T as identified from Figure 4 (in Figure 6, ρz(x)

Figure 6. ρz(x) at ρ = 0.732 and various T values in the range of stability of the two Cassie−Baxter states. The tangent angle α is also reported. Note that the peculiar form of the meniscus at T = 0.76 is due to the existence of two specular asymmetric configurations, whose average is represented in the panel. A similar phenomenon is observed more evidently at higher z, as discussed in detail in Mechanism of the CB/W Transition.

are shown for the nominal density ρ = 0.732). The results are reported in Figure 7. They show that α has a discontinuity when changing phase from the CB1 to the CB2 state. Moreover, the dependence of α on ρ and T when the system is in the CB1 state is negligible: its value is always ∼90°. On the

Figure 7. α as a function of the density and temperature. White regions on the graph correspond to ρ−T values for which we did not perform any simulation. 10768

dx.doi.org/10.1021/la3018453 | Langmuir 2012, 28, 10764−10772

Langmuir

Article

states, and Figure 9B for the same field along the transition path). This observation suggests that the atomistic mechanism may follow unexpected asymmetric paths. To make this argument independent of subjective observation, we apply cluster analysis59 to the instantaneous density fields generated by RMD-PT in order to group similar configurations into clusters (an extended description of this approach is given in Appendix 1). A transition path is then extracted as the sequence of clusters at successive values of z that are the closest to each other. The concepts of “similarity” and “proximity” need the introduction of the notion of distance between two instantaneous density fields and between two clusters. The reader interested in the algorithmic details is referred to the Appendix 1; for the moment, we only mention that we measured the distance between two configurations indexed by k and k′, corresponding in general to two different z, in terms of the L2-norm of the density field, Δρ̃kk' zz′ = [∫ dx (ρ̃kz(x) − ρ̃k'z′(x))2]1/2. Once clusters are formed at each z, a representative of each cluster is selected, the medoid, that is the element of minimal average dissimilarity with all the other elements of the cluster. A ρ̃z(x) vs z trajectory is built starting from one medoid at z = 1800, that corresponds to the transition state (maximum of the free energy F(z)). For this medoid, we identify the closest medoid at the next z value. This procedure is iterated for any pair of contiguous z, moving toward CB2 and toward W until the complete sequence of clusters connecting the two states is constructed. The entire process is repeated for another cluster at the transition state, until the paths starting from all these clusters have been identified. We applied this analysis to the simulations at ρ = 0.732 and T = 0.81. As a first remark, we mention that we were able to identify two well-separated clusters for z ∈ [1600, 1800], while at lower (z ∈ [1200, 1600)) and higher (z ∈ (1800, 2000)) z, our data are better represented by a single cluster. This confirms that the transition can proceed along two different paths. These paths initially coincide and then split at intermediate z, merging back at higher values (see Figure 9A). For comparison, in Figure 9B we report the ρz(x) vs z path, i.e., the one obtained by averaging over all the states sampled at a given z values. The instantaneous density fields ρ̃z(x) show that initially (1200 ≤ z ≤ 1500) the meniscus descends flatly, then (1600 ≤ z ≤ 1800) it bends forming a vapor bubble in one corner of the groove (symmetry breaking), and finally (z ≥ 1900) the bubble is absorbed and the groove filled. The transition state, occurring at z = 1800 at this ρ−T, corresponds to configurations in which a bubble is formed in one of the corners. This mechanism is significantly different from the flat meniscus considered in continuum literature,21,40 that compares well with Figure 9 only at low z. Another mechanism has also been proposed, namely, the “sag” transition,41 where the triple line remains pinned at the edges of the groove and the liquid/vapor interface bends symmetrically toward the bottom. In the present simulations a similar symmetric shape is obtained by averaging over the two specular atomistic paths (Figure 9B) at z ≈ 1800. A possible macroscopic explanation of the asymmetric transition can be obtained considering the 2D problem consisting of a square box partially filled with a liquid in equilibrium with its vapor (see also Supporting Information). This problem is similar to the one considered in this work if one assumes that there is translational invariance along the z direction (see Figure 2). If the liquid has a Young contact angle

Figure 8. Free energy barrier ΔF† of the CB2 → W and W → CB2 transition as a function of T for ρ = 0.725 and 0.732. ΔF† is reported in kBT units. The dashed line in the two panels denotes the estimated CB/W transition temperature.

our simulations, as we do not have trajectories going from CB2 to W. In fact, RMD-PT provides sets of configurations extracted from the probability density function conditional to ε(r) = z at fixed z (see Methods). These configurations represent highly probable states at a fixed value of z. We note that, for some z, there exist, at least from visual inspection, two classes of specular atomistic configurations that break the reflectional symmetry with respect to the central y−z plane. An example is provided in Figure 9A where the plots refer to

Figure 9. CB2 to W transition paths at T = 0.81 and ρ = 0.732. The z values range from 1200 to 1900. In panel (A) the smoothed instantaneous field ρ̃z(x) is reported along the two specular paths described in the text. Panel (B) provides the ensemble averaged field ρz(x) at the corresponding z.

smoothed instantaneous density fields ρ̃z(x).52 A transition between the two asymmetric states was never observed in a run without PT. This is an indication that they represent free energy minima of unknown collective variables (belonging to the subspace orthogonal to ε(r)) that are separated by large free energy barriers. Only using PT, these free energy barriers can be overcome, and the symmetry of the ensemble averaged ρz(x) is recovered (see Figures 5 and 6 for ρz(x) in equilibrium 10769

dx.doi.org/10.1021/la3018453 | Langmuir 2012, 28, 10764−10772

Langmuir

Article

θY = 90° (i.e., when the solid/vapor and solid/liquid interface energies are the same), and under the usual assumptions (e.g., that the surface energy and the chemical potential is constant at all levels of filling of the box), the minimum of the free energy is attained by minimizing the liquid/vapor interface (a line in the 2D case). At low filling levels, the minimizing configuration is a straight line (see Figure 10A), while for higher liquid

the transition to the W state will happen spontaneously, as the metastable CB state does not exist any longer. From a simulative point of view, metastabilities confirm that the CB-W transition is indeed a rare event. This implies that the study of nanoscale wetting by atomistic simulations requires sophisticated sampling methods. Concerning the mechanism of the transition, we found that it follows an asymmetric path, with the liquid forming a finger on one side of the groove and the vapor forming a bubble on the opposite side, a fact that could be qualitatively interpreted in terms of a simple macroscopic argument based on the minimization of surface energy.

■

APPENDIX 1: CLUSTER ANALYSIS FOR THE IDENTIFICATION OF POSSIBLE TRANSITION PATHS In this section, we describe in more detail the calculation of the CB-to-W transition path. The first point to stress is that the instantaneous field ρ̂z(x) is noisy and this could make the convergence of the clustering algorithm more difficult and its results unreliable. To solve this problem, we replace ρ̂z(x) with a smooth equivalent obtained by averaging over a time window centered around the time corresponding to the relevant 52,60 configuration. In practice, ρ̃z(x) = 1/τ∫ t+τ/2 t−τ/2 ds ρ̂ z(x, r(s)). Using this smooth field, a classification of the instantaneous densities is performed with the aid of clustering analysis.59 This class of techniques is capable of performing a rapid classification of large data sets and to provide quantitative information about the separation between data and clusters. The principle behind any clustering algorithm is associating different data, based on a conveniently defined distance. For two configurations indexed by i and j, we use hereafter the distance based on the L2-norm of the density field Δρ̃ijz = [∫ dx(ρ̃iz(x) − ρ̃jz(x))]1/2. In practice, the observable ρ̃z(x) is defined on a discretization of the simulation box, and the integral is substituted by the corresponding sum over the discretized cells. Local density values lower than 0.1 are filtered out, in order to eliminate the noise due to the vapor phase. We performed all cluster analysis using the R environment.61 The first insight into possible data classification comes from hierarchical clustering,62 that, starting from the single data, at each iteration merges the closest two clusters. The proximity between clusters is evaluated with the Ward’s method.59 Data are organized in a dendrogram, where the root of the tree is the all-comprehensive cluster and the leaves represent the single data. The length of the tree branches measures the proximity between the merged clusters. The picture emerging from this first classification is that data at 1600 ≤ z ≤ 1800 are likely to be organized in two clusters, while for lower and higher z, it is difficult to distinguish among clusters. Informed by hierarchical clustering, we employ the PAM (partitioning around medoids) algorithm63 which attempts to form a user-selected number of clusters, and yields a measure of the goodness of a clustering (the silhouette coef ficient going from 0 to 1). In addition, PAM selects a representative of each cluster, the medoid, chosen as the closest element to the cluster center. By requiring formation of two clusters, we obtain a silhouette coefficient higher than 0.5 for 1600 ≤ z ≤ 1800, indicating a robust structure of the clusters (this physically corresponds to fluid wetting either the left or the right corner). For the other z, the clustering structure is weak, with a silhouette coefficient ∼0.25. We form the succession of medoids representing the transition, as in Figure 9A, by evaluating the closest medoids for contiguous values of z. We

Figure 10. Sketch of a 2D groove partly filled with a liquid in equilibrium with a vapor: (A) flat liquid/vapor interface, (B) quarter of circle liquid/vapor interface. For R > 2l/π, the flat interface has the lowest free energy; for R < 2l/π, the most stable interface is (B).

volumes in the box, the minimal free energy corresponds to a quarter of circle occupying one corner (Figure 10B). The transition from the flat to the quarter of circle case occurs when the radius of the latter is lower than 2l/π, where l is the width of the box (see Figure 2 for comparison). Summarizing, when θ ≈ 90°, the flat-to-curved meniscus transition, observed in the data shown in Figure 9 at high values of z, may be explained as a minimization of the liquid/vapor interface area. This shows that, at least at a qualitative level, macroscopic theories are able to describe the CB/W transition occurring on nanopatterned surfaces once unneeded assumptions, such as the symmetry of the density field and the pinning of the triple line, are lifted.

■

CONCLUSIONS In this paper, we studied the phase diagram with respect to temperature and (nominal) density of a simple (LJ) bulk liquid contained between two walls featuring one groove. Given the size of the groove in terms of atomic radius, our geometry corresponds to a nanocorrugation. For such a system, the simple two-state (Cassie−Baxter and Wenzel) representation is insufficient as a second Cassie−Baxter state characterized by a large meniscus curvature has been identified. We used the angle formed by the tangent of the meniscus and the wall to characterize the system dependency on temperature and density, finding that higher temperatures tend to flatten the meniscus down, while the density has the opposite effect. This angle is of great importance as our simulations show that a flat meniscus must be reached before the Cassie−Baxter/Wenzel transition takes place, in accordance with Gibbs’ criterion. Cassie−Baxter/Wenzel metastability has also been investigated. The free energy barrier separating these two states was found to largely exceed the thermal energy, confirming that liquids can be trapped in these metastable states. The Cassie− Baxter state and CB/W metastabilities are supported even by the slightly hydrophilic surface considered (θY ≲ 90°). This evidence suggest that surfaces having the simple geometry considered here may be realized (for instance, by photolithographic techniques on silicon wafers) to support the superhydrophobic CB state over a wide range of temperatures and densities (or, equivalently, pressures) and for a variety of liquids, provided that the Young contact angle is approximately equal to or larger than 90°. In addition, we provide evidence that there exist limit temperatures and densities beyond which 10770

dx.doi.org/10.1021/la3018453 | Langmuir 2012, 28, 10764−10772

Langmuir

Article

(4) Johnson, R. E.; Dettre, R. H. In Contact Angle, Wettability, and Adhesion, Fowkes, F. M., Ed.; Advances in Chemistry Series; ACS Publications: Washington, DC, 1964; Vol. 43, Chapter 7, pp 112−135. (5) Dettre, R. H.; Johnson, R. E. In Contact Angle, Wettability, and Adhesion, Fowkes, F. M., Ed.; Advances in Chemistry Series; ACS Publications: Washington, DC, 1964; Vol. 43, Chapter 8, pp 136−144. (6) Joanny, J.; De Gennes, P. A Model for Contact Angle Hysteresis. J. Chem. Phys. 1984, 81, 552−562. (7) Koishi, T.; Yasuoka, K.; Fujikawa, S.; Zeng, X. Measurement of Contact-Angle Hysteresis for Droplets on Nanopillared Surface and in the Cassie and Wenzel States: A Molecular Dynamics Simulation Study. ACS Nano 2011, 5, 6834−6842. (8) Cassie, A. B. D.; Baxter, S. Wettability of Porous Surfaces. Trans. Faraday Soc. 1944, 40, 546−551. (9) In their derivation, Cassie and Baxter considered coplanar liquid− vapor and solid−liquid interfaces and used this configuration to estimate the effective contact angle θCB; see Figure 1A. In most real cases, however, the form of the liquid−vapor interface is much more complicated and θCB needs a correction.10 Since we are not interested directly in effective contact angles, we will refer to the CB state somewhat loosely, and include in its definition all (meta)stable states where vapor is trapped inside the groove. (10) Milne, A. J. B.; Amirfazli, A. The Cassie Equation: How It is Meant to be Used. Adv. Colloid Interface Sci. 2012, 170, 48−55. (11) Wenzel, R. N. Resistance of Solid Surfaces to Wetting by Water. Ind. Eng. Chem. 1936, 28, 988−994. (12) In a fluid experiencing a mean flow, the slip length Ls is defined as the distance from the liquid/solid interface at which, in the linear regime approximation, the velocity field of the liquid gets to zero. The convention is that the slip length is positive when this condition is met in the wall side, and negative otherwise. Therefore, to higher values of Ls correspond higher velocities of the fluid at the wall. (13) Cottin-Bizonne, C.; Cross, B.; Steinberger, A.; Charlaix, E. Boundary Slip on Smooth Hydrophobic Surfaces: Intrinsic Effects and Possible Artifacts. Phys. Rev. Lett. 2005, 94, 056102. (14) Chinappi, M.; Casciola, C. M. Intrinsic Slip on Hydrophobic Self-Assembled Monolayer Coatings. Phys. Fluids 2010, 22, 042003. (15) Lee, C.; Choi, C.-H.; Kim, C.-J. Structured Surfaces for a Giant Liquid Slip. Phys. Rev. Lett. 2008, 101, 064501. (16) Bahners, T.; Textor, T.; Opwis, K.; Schollmeyer, E. Recent Approaches to Highly Hydrophobic Textile Surfaces. J. Adhes. Sci. Technol. 2008, 22, 285−309. (17) Liu, Y.; Chen, X.; Xin, J. H. Can Superhydrophobic Surfaces Repel Hot Water? J. Mater. Chem. 2009, 19, 5602−5611. (18) Rothstein, J. P. Slip on Superhydrophobic Surfaces. Annu. Rev. Fluid Mech. 2010, 42, 89−109. (19) Zhang, X.; Shi, F.; Niu, J.; Jiang, Y.; Wang, Z. Superhydrophobic Surfaces: From Structural Control to Functional Application. J. Mater. Chem. 2008, 18, 621−633. (20) Extrand, C. W. Criteria for Ultralyophobic Surfaces. Langmuir 2004, 20, 5013−5018. (21) Patankar, N. A. Transition between Superhydrophobic States on Rough Surfaces. Langmuir 2004, 20, 7097−7102. (22) Bhushan, B.; Nosonovsky, M.; Jung, Y. C. Towards Optimization of Patterned Superhydrophobic Surfaces. J. R. Soc. Interface 2007, 4, 643−648. (23) Lee, J. B.; Gwon, H. R.; Lee, S. H.; Cho, M. Wetting Transition Characteristics on Microstructured Hydrophobic Surfaces. Mater. Trans. 2010, 51, 1709−1711. (24) Afferrante, L.; Carbone, G. Microstructured Superhydrorepellent Surfaces: Effect of Drop Pressure on Fakir-State Stability and Apparent Contact Angles. J. Phys.: Condens. Matter 2010, 22, 325107. (25) Koishi, T.; Yasuoka, K.; Fujikawa, S.; Ebisuzaki, T.; Zeng, X. Coexistence and Transition between Cassie and Wenzel State on Pillared Hydrophobic Surface. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 8435−8440. (26) Savoy, E. S.; Escobedo, F. A. Molecular Simulations of Wetting of a Rough Surface by an Oily Fluid: Effect of Topology, Chemistry

start from z = 1800, where a well-defined clustering exists (and incidentally corresponds to the transition state), and apply the algorithm in the CB2 and W directions. We notice that for 1200 ≤ z < 1600 and z ≥ 1900 the separation between clusters at the same z decreases. This is consistent with the weak clustering of the symmetric fields at low/high z, which actually form a single cluster. The transition path analysis we use in this investigation is based on the assumption that the transition process in the discretized z-space can be described in terms of a discrete-state stochastic dynamics, with an associated stochastic matrix whose only nonzero elements are those corresponding to states belonging to previous/successive z values. Moreover, the nonzero elements of the stochastic matrix quickly decrease with the corresponding Δρ̃ijzz′, which is computed over the medoids of the ith and jth clusters at z and z′, respectively. A possible justification for these hypotheses is the following. The stochastic matrix elements between clusters with same z are small as these clusters are separated by large free energy barriers. This is, indeed, confirmed by our simulations. In fact, due to the presence of these barriers, only using PT we were able to correctly sample the conditional invariant distribution at a given z value. As for the dependency of the stochastic matrix on the norm of the density difference Δρ̃ijzz′, this is explained by considering that ρ̃iz(x) at a given z is assumed to be the evolution of ρ̃iz′(x) at the previous z value via a continuous dynamics. Thus, we expect that the two fields are not too dissimilar and, in turn, the probability to have a path connecting two states with a large density difference is low.

■

ASSOCIATED CONTENT

S Supporting Information *

Simulation Details. Contact angles of a LJ liquid on a FCC flat surface. Minimal Surface Energy inside a Rectangular Groove. This material is available free of charge via the Internet at http://pubs.acs.org/.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors wish to acknowledge L. Ferraro for his support in setting up the computational environment for the present simulations. The authors wish to acknowledge the Consorzio per le Applicazioni del Supercalcolo Per Università e Ricerca (CASPUR) and SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational resources. This research was supported by a Marie Curie Intra European Fellowship within the seventh European Community Framework Programme.

■

REFERENCES

(1) Gao, L.; McCarthy, T. J. A Perfectly Hydrophobic Surface (θA/θR = 180°/180°). J. Am. Chem. Soc. 2006, 128, 9052−9053. (2) Quéré, D. Wetting and Roughness. Annu. Rev. Mater. Res. 2008, 38, 71−99. (3) Reyssat, M.; Pépin, A.; Marty, F.; Chen, Y.; Quéré, D. Bouncing Transitions on Microtextured Materials. Europhys. Lett. 2006, 74, 306−312. 10771

dx.doi.org/10.1021/la3018453 | Langmuir 2012, 28, 10764−10772

Langmuir

Article

each ρ−T point. Then, at each z, we computed the error on the F(z) by the block average technique. This approach produces a more accurate estimate of the error as compared to the error propagation method, often used in the literature. In fact, the error propagation method systematically overestimates it, giving rise to an unrealistic quick divergence along the curve. (50) The “nominal” density is defined as ρ = N/V, where N is the number of liquid particles and V is the volume available to the liquid. We fix the density rather than the pressure as this allows to significantly reduce the computational cost of our simulations. In the thermodynamic limit, the equivalence of the ensembles guarantees that the fixed pressure and fixed density results are equivalent. However, the finite number of particles in our computational sample might lead to non negligible differences. We plan to investigate the possible effects of pressure in a future work. (51) Hansen, J.-P.; Verlet, L. Phase Transitions of the Lennard-Jones System. Phys. Rev. 1969, 184, 151−161. (52) Following Kirkwood,64 the density field is defined as ρ(x) = ⟨ρ̂(x,r )⟩ = ⟨∑i=1,N δ(ri − x)⟩, where x denotes a point in the ordinary 9 3 space, the sum runs over the atoms in the sample and ⟨·⟩ stands for an ensemble average. The symbol ρ̂z(x, r) denotes an instantaneous density field at a given z, while ρz(x) = ⟨ρ̂(x, r)⟩z denotes the expectation value of the density operator over the ensemble at ε(r) = z, i.e., over the conditional probability density function w(r)δ(ε(r)− z)/∫ dr w(r)δ(ε(r)−z). This conditional expectation value can be obtained by restrained MD.65 The instantaneous density field ρ̂z(x, r) is noisy and its smoothed version ρ̃z(x,r) = 1/τ∫ t+τ/2 t−τ/2 ds ρ̂z(x,r(s)) is used throughout. (53) Gibbs, J. W. The Scientific Papers of J. Willard Gibbs, Vol. 1: Thermodynamics; Dover Publications: New York, 1961; p 434. (54) Oliver, J.; Huh, C.; Mason, S. Resistance to Spreading of Liquids by Sharp Edges. J. Colloid Interface Sci. 1977, 59, 568−581. (55) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Transition Path Sampling: Throwing Ropes over Rough Mountain Passes, in the Dark. Annu. Rev. Phys. Chem. 2002, 53, 291−318. (56) Marmur, A. Soft Contact: Measurement and Interpretation of Contact Angles. Soft Matter 2006, 2, 12−17. (57) Boreyko, J.; Chen, C.-H. Restoring Superhydrophobicity of Lotus Leaves with Vibration-Induced Dewetting. Phys. Rev. Lett. 2009, 103, 174502. (58) Dorrer, C.; Rühe, J. Some Thoughts on Superhydrophobic Wetting. Soft Matter 2009, 5, 51−61. (59) Tan, P.-N.; Steinbach, M.; Kumar, V. Introduction to Data Mining; Addison-Wesley, 2006. (60) Care must be paid to this window time average to avoid to include configurations coming from different PT trajectories after a successful swapping. In fact, this will not correspond to a smoothing but to an (insufficient) ensemble average. (61) R Development Core Team; R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing; Vienna, Austria, 2010. (62) Murtagh, F. COMPSTAT Lectures 4; Physica-Verlag, 1985. (63) Kaufman, L.; Rousseeuw, P. Finding Groups in Data: An Introduction to Cluster Analysis; J. Wiley and Sons, Inc., 1990. (64) Irving, J. H.; Kirkwood, J. G. The Statistical Mechanical Theory of Transport Processes. IV. The Equations of Hydrodynamics. J. Chem. Phys. 1950, 18, 817. (65) Orlandini, S.; Meloni, S.; Ciccotti, G. Hydrodynamics from Statistical Mechanics: Combined Dynamical-NEMD and Conditional Sampling to Relax an Interface between Two Immiscible Liquids. Phys. Chem. Chem. Phys. 2011, 13, 13177.

and Droplet Size on Wetting Transition Rates. Langmuir 2012, 28, 3412−3419. (27) Cottin-Bizonne, C.; Barrat, J.-L.; Bocquet, L.; Charlaix, E. LowFriction Flows of Liquid at Nanopatterned Interfaces. Nat. Mater. 2003, 2, 237−240. (28) Lafuma, A.; Quéré, D. Superhydrophobic States. Nat. Mater. 2003, 2, 457−460. (29) Martines, E.; Seunarine, K.; Morgan, H.; Gadegaard, N.; Wilkinson, C. D. W.; Riehle, M. O. Superhydrophobicity and Superhydrophilicity of Regular Nanopatterns. Nano Lett. 2005, 5, 2097−2103. (30) Here, we refer to “equilibrium phase diagram” to distinguish it from the one that can be obtained from experiments on nonequilibrium systems phase diagram. (31) Orlandini, S.; Meloni, S.; Ciccotti, G. Combining Rare Events Techniques: Phase Change in Si Nanoparticles. J. Stat. Phys. 2011, 145, 812−830. (32) Zisman, W. In Contact Angle, Wettability, and Adhesion Fowkes, F. M., Ed.; Advances in Chemistry Series; ACS Publications: Washington, DC, 1964; Vol. 43, Chapter 1, pp 1−51. (33) Tuteja, A.; Choi, W.; Ma, M.; Mabry, J. M.; Mazzella, S. A.; Rutledge, G. C.; McKinley, G. H.; Cohen, R. E. Designing Superoleophobic Surfaces. Science 2007, 318, 1618−1622. (34) Reyssat, M.; Yeomans, J. M.; Quéré, D. Impalement of Fakir Drops. Europhys. Lett. 2008, 81, 26006. (35) Kusumaatmaja, H.; Blow, M. L.; Dupuis, A.; Yeomans, J. M. The Collapse Transition on Superhydrophobic Surfaces. Europhys. Lett. 2008, 81, 36003. (36) Vrancken, R. J.; Kusumaatmaja, H.; Hermans, K.; Prenen, A. M.; Pierre-Louis, O.; Bastiaansen, C. W. M.; Broer, D. J. Fully Reversible Transition from Wenzel to Cassie-Baxter States on Corrugated Superhydrophobic Surfaces. Langmuir 2010, 26, 3335−3341. (37) Whyman, G.; Bormashenko, E. How to Make the Cassie Wetting State Stable? Langmuir 2011, 27, 8171−8176. (38) Chandler, D. Interfaces and the Driving Force of Hydrophobic Assembly. Nature 2005, 437, 640−647. (39) Alberti, G.; De Simone, A. Wetting of Rough Surfaces: A Homogenization Approach. Proc. R. Soc. London, Ser. A 2005, 461, 79− 97. (40) Ishino, C.; Okumura, K.; Quéré, D. Wetting Transitions on Rough Surfaces. Europhys. Lett. 2004, 68, 419. (41) Patankar, N. A. Consolidation of Hydrophobic Transition Criteria by Using an Approximate Energy Minimization Approach. Langmuir 2010, 26, 8941−8945. (42) Maragliano, L.; Vanden-Eijnden, E. A Temperature Accelerated Method for Sampling Free Energy and Determining Reaction Pathways in Rare Events Simulations. Chem. Phys. Lett. 2006, 426, 168−175. (43) Swendsen, R. H.; Wang, J. S. Replica Monte Carlo Simulation of Spin-Glasses. Phys. Rev. Lett. 1986, 57, 2607−2609. (44) Earl, D. J.; Deem, M. W. Parallel Tempering: Theory, Applications, and New Perspectives. Phys. Chem. Chem. Phys. 2005, 7, 3910−3916. (45) Kirkwood, J. C. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3, 300−313. (46) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (47) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R.; Parrinello, M. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972. (48) Meloni, S.; Rosati, M.; Colombo, L. Efficient Particle Labeling in Atomistic Simulations. J. Chem. Phys. 2007, 126, 121102. (49) We computed the error bars by applying the block average method to the results of the thermodynamic integration. First, we divided the samples at fixed z in smaller sets and computed the free energy gradient over these sets. From these sets of mean forces, by thermodynamic integration, we obtained several free energy curves at 10772

dx.doi.org/10.1021/la3018453 | Langmuir 2012, 28, 10764−10772