Sara
Marchio
a,
Simone
Meloni
*ab,
Alberto
Giacomello
*a and
Carlo Massimo
Casciola
a
aDipartimento di Ingegneria Meccanica e Aerospaziale, Università di Roma Sapienza, Via Eudossiana 18, 00184 Roma, Italy. E-mail: alberto.giacomello@uniroma1.it
bDipartimento di Scienze Chimiche e Farmaceutiche (DipSCF), Universitá degli Studi di Ferrara (Unife), Via Luigi Borsari 46, I-44121, Ferrara, Italy. E-mail: simone.meloni@unife.it
First published on 4th November 2019
Hydrophobic (nano)textured surfaces, also known as superhydrophobic surfaces, have a wide range of technological applications, including in the self-cleaning, anti-moisture, anti-icing, anti-fogging and friction/drag reduction fields, and many more. The accidental complete wetting of surface textures, which destroys superhydrophobicity, and the opposite process of recovery are two crucial processes that can prevent or enable the technological applications mentioned before. Understanding these processes is key to designing surfaces with tailored wetting and recovery properties. However, recent experiments have suggested that the currently available theories are insufficient for describing the observed phenomena. In this work we offer a dynamic picture of these processes beyond the state of the art showing that the key ingredient determining the experimental behavior is the inertia of the liquid in the wetting and dewetting processes, which is neglected in microscopic and macroscopic quasi-static theories inspired by the classical nucleation theory. The present findings are also important for other related phenomena, such as heterogeneous cavitation, where vapor/gas bubbles form at surface asperities, condensation, dynamics of the triple line, micelle formation, etc.
Textured surfaces can be fabricated using several top down and bottom up techniques (see ref. 14 for a comprehensive review of fabrication techniques). Among others, a widely used technique is lithography in its different forms: photolithography, soft lithography, nanoimprint lithography, electron beam lithography, X-ray lithography, and colloidal lithography. In lithography, textured surfaces are prepared by copying the information from a master on a substrate using a mask. Very promising is thermal scanning probe lithography,15 enabling the fabrication of surfaces with sub-10 nm features. Another emerging approach is self-assembly, in which the formation of structures, aggregates and films is a spontaneous process driven by molecular and/or supramolecular forces. Recently, self-assembly of block-copolymers has been used to fabricate superhydrophobic surfaces with ∼10 nm textures.16
In addition to the Cassie–Baxter state a second state exists, known as the Wenzel state,17 in which the liquid completely wets the corrugations of the hydrophobic surface. In the Wenzel state the liquid/solid contact area is much larger and the superhydrophobic properties are lost. The Cassie–Baxter/Wenzel transition (wetting – Fig. 1) can be induced by changes of pressure and/or temperature.18,19 The fragility of the Cassie–Baxter state and the difficulty of recovery (Wenzel/Cassie–Baxter transition20–22) have hindered the use of superhydrophobic surfaces in practical applications. Thus, intense experimental and theoretical research has been devoted to the investigation of the wetting and recovery mechanisms.20,22–38
In their seminal work, Lafuma and Quéré22 have investigated the variation of the contact angle hysteresis as a function of the pressure of a liquid droplet deposited between two textured hydrophobic surfaces. A large and sudden increase of the contact angle hysteresis is a measure of the transition from the Cassie–Baxter to the Wenzel state. In particular, Lafuma and Quéré have shown that upon relaxing the pressure the contact angle hysteresis does not recover the original value; that is, the transition is irreversible.
The experimental investigation of the wetting mechanism in submerged surfaces poses challenges for following the evolution of the state of the system with the change of external conditions. Optical techniques have been employed for this purpose, such as light reflection stereo microscopy (see, e.g., ref. 39) and transmission diffraction (see, e.g., ref. 20, 40 and 41). In contrast with these indirect techniques, confocal microscopy directly follows the evolution of the liquid–air meniscus,24,42,43 which makes it an optimal tool for investigating the wetting and recovery mechanism. In particular, Duan and coworkers used confocal microscopy to follow the evolution of the water meniscus while the pressurized liquid enters into the cylindrical pores decorating a hydrophobic surface.24,43 They have shown that water intrudes the pores with a symmetric, almost flat, meniscus. Indeed, also an asymmetric wetting path has been observed in these experiments, which consists in the formation of a gas bubble at the bottom of corrugations which is then absorbed. The observation of this second wetting path has been attributed to the presence of impurities accumulating on the side walls of corrugations that pin the meniscus during its advancement.
Under the conditions of interest for experiments and technological applications the wetting and recovery transitions are often characterized by (relatively) large free-energy barriers (ΔΩ†w and ΔΩ†r, respectively) separating the initial and final states. In the presence of barriers larger than the thermal energy kBT (kB Boltzmann constant and T temperature) the transition time scales exponentially with ΔΩ†w/r:44,45
τw/r = τ0w/r exp(ΔΩ†w/r/kBT). | (1) |
Because of this long transition time, largely exceeding the timescale accessible to brute force (continuum and/or atomistic) simulations, special techniques are necessary to investigate the wetting and recovery transitions. Apart from notable exceptions,46,47 the wetting and recovery of textured surfaces has been studied via quasi-static methods, such as umbrella sampling,31 restrained molecular dynamics,27,37,48 string and nudged elastic band methods.32,33,36,48–51 These methods assume that the process is slow and the system is at the local equilibrium all along wetting or recovery. In other words, these methods neglect dynamic effects, such as inertia. For a set of textured surfaces, e.g. surfaces with circular and square pores, ridges and pillars, the most probable wetting path predicted by quasi-static approaches – the one associated with the lowest free-energy barrier – breaks the symmetry of the system; for example, it is characterized by the formation of a gas bubble in the corner of the corrugations.27–29,48,52 The comparison with complementary quasi-static macroscopic theories and simulations28,32 confirms that the characteristic length of surface corrugations does not alter the overall picture of the wetting mechanism, which remains asymmetric. These results confirm that quasi-static theoretical predictions are at odds with the recent experiments discussed above.
The objective of this work is to reconcile the mismatch between theory and experiments, i.e. to identify the neglected terms that are responsible for the symmetric wetting and to extend existing theories to predict a wetting path consistent with experimental observations. This extended framework is expected to be relevant also for related phenomena, such as homogeneous and heterogeneous bubble nucleation, which hardly conforms to a quasi-static picture.53,54 To achieve our aim we perform molecular dynamics (MD) simulations that, as discussed in more detail in the next section, embodies all relevant phenomena affecting the evolution of the solid–liquid–gas three phase system. The use of atomistic simulations, however, imposes some constraint on the size of textures that can be simulated, between 1 and 10 nm. Nevertheless, as discussed in more detail in the next section, fluid mechanics arguments based on dimensionless numbers, namely the Womersley number, confirm that our conclusions can be extended to larger textures of experimental and applicative interest.
We investigate a wide range of thermodynamic conditions. In particular, we run simulations at pressures corresponding to high and low wetting barriers associated, respectively, with long and short transition times. The former corresponds to ambient conditions, of practical interest for most of the common textured hydrophobic surfaces, and the latter to more extreme conditions used in wetting experiments.
The present results are of interest from both the fundamental and the application points of view: in addition to revealing paramount dynamic effects in nucleation, such a theory opens a way to design surfaces with tailored properties that can better resist wetting32,55,56 and enable facile recovery.31,52 Similar dynamic effects are expected to play a role in the broader context of heterogeneous nucleation, liquid fronts moving on actual surfaces, self-assembly, micelle formation, etc. Moreover, wetting and recovery, liquid intrusion in and extrusion from lyophobic pores, is important for the emerging field of mechanical energy storage and energy dissipation using nanometric hydrophobic porous systems, such as mesoporous silica, e.g. MCM-4157 (pore size ∼2–6.5 nm), or metal–organic frameworks (MOFs), e.g. ZIF-858,59 (pore size ∼1.2 nm).
We remark that our simulations allow reconciling the experimental and theoretical pictures of wetting: when the free-energy barrier is sufficiently low – close to the values at which the experimental transition is expected to occur – dynamic effects dominate and the meniscus advances inertially in the pores, preserving the initial symmetric shape. In addition, our results support the interpretation that the asymmetric wetting observed from time to time in experiments is due to the deposition of impurities on the textured walls. Finally, in contrast with the quasi-static picture,28,31–33,36,37,48,50,51 the present results show that dynamic effects are strongly dependent on the liquid pressure.
In the present work we use RMD to provide (quasi-static) reference results for comparison. Here, we control (restrain) the number of fluid particles in the pore N (Fig. 2A). A continuum fluid dynamics analysis (Navier–Stokes equations60) may help understanding the limitations of quasi-static approaches. In atomistic simulations the wetting/recovery liquid velocity can be recast into the form u = mṄ/(Aρ), where Ṅ is liquid particles’ wetting/recovery rate, A is the area of the pore, and m and ρ are the particle mass and the liquid density, respectively. In quasi-static approaches Ṅ is set to 0 by the restraint at any point along the process and so is the liquid velocity u. This amounts to neglecting inertial over viscous forces, as quantified by the Womersley number
α2 = 2π Sr Re = ρL2/μτ = mṄ/μL = 0, | (2) |
Our FFS simulations show that, when wetting and recovery take place on the experimental timescale, α2 reaches values up to 0.33 (see below), suggesting that inertia may play an important role in modeling these processes. We remark that, when dimensionless numbers are invaluable to quantify inertial effects, (mean-field) Navier–Stokes equations cannot be used to simulate thermally activated processes, where thermal fluctuations are the key to overcome free-energy barriers. The atomistic simulations discussed in the following incorporate such thermal fluctuations and enabled us to assess the relevance of inertia in the wetting/recovery process under varying thermodynamic conditions – specifically liquid pressure.
We first consider a Lennard–Jones (LJ) liquid wetting a surface made of LJ particles and featuring a 10 × 10 ridge (∼3.5 nm × 3.5 nm if the liquid were water), i.e. a square cavity extending infinitely in the direction orthogonal to the plane of Fig. 1. This geometry has been extensively investigated in previous theoretical28,38,61,62 and computational27,30,48 studies because despite its simplicity it embodies all the ingredients governing the wetting/recovery transition.61 The fluid and solid particles interact via the modified LJ potential
(3) |
We performed RMD simulations, i.e. without inertia, of the wetting and recovery processes at different liquid pressures, corresponding to different values of the free energy difference ΔΩCB/W between the Cassie–Baxter and Wenzel states and of the wetting/recovery barrier ΔΩ†w/r: higher pressures favor the stability of the Wenzel state and reduce the wetting barrier and vice versa (see Fig. 3 and Table 1). In particular, at P = −0.168 and P = −0.08 Cassie–Baxter is the stable state and Wenzel is the metastable one; the situation is reversed at P = −0.005, 0.01, and 0.035: Wenzel is the stable state and Cassie–Baxter is the metastable one. According to eqn (1), only when the barrier is of the order of few kBT does the transition take place on the experimental timescale (milliseconds to seconds). Thus, considering the wetting and recovery barriers reported in the table, our RMD results suggest that, among the pressures considered, wetting readily takes place at P = 0.035 and recovery at P = −0.168.
P | ΔΩCB/W | ΔΩ†w | ΔΩ†r |
---|---|---|---|
-0.168 (−3.8) | −130 | 133 | 3 |
-0.08 (−1.8) | -74 | 91 | 17 |
-0.005 (−0.1) | 21 | 17 | 38 |
0.01 (0.2) | 31 | 12 | 43 |
0.035 (0.8) | 57 | 0 | 57 |
Consistent with literature results,27,28,32 the RMD wetting transition begins with liquid depinning from the top corners and entering in the cavity with an (almost) flat meniscus (see Fig. 4A); in correspondence of the transition state (the maximum of the free-energy profile, N ∼ 450) the meniscus bends to form a vapor bubble in one of the two bottom corners of the cavity. Then, the bubble shrinks until it disappears when the system reaches the Wenzel state.
Fig. 4 (A) Sequence of snapshots along the wetting (left-to-right) and recovery (right-to-left) paths as obtained by restrained molecular dynamics.27,37,69 During the wetting the liquid initially enters in the pore with a flat meniscus, and then forms a bubble in a corner and finally the bubble is absorbed and the meniscus touches the bottom wall. The quasi-static process is reversible and thus the recovery path is the reverse of the wetting path. (B) Pairs of (ΔN, N) values sampled in RMD simulations of the wetting/recovery process at P = −0.08. Analogous graphs for the other pressures are shown in Fig. SM2.† The black dashed line, introduced to help the reader to follow the wetting/recovery path, is obtained by connecting the average value of ΔN at each N, 〈ΔN〉N; in the region in which the graph is split into two parts (N > 450) 〈ΔN〉N is computed within the domain of each mode. Consistent with the snapshots of panel (A), in the early part of the wetting ΔN values are concentrated about 0. At N ∼ 450 one observes a sharp change, with sizably negative and positive ΔN values. |
A more quantitative determination of the asymmetry of the liquid front along the wetting process is provided by ΔN = N1 − N2vs. N, where N1 and N2 are the numbers of liquid particles in the left and right halves of the cavity, respectively (Fig. 2B). When the meniscus is flat ΔN ∼ 0, and when a bubble is formed in a corner ΔN is either sizably negative or positive. The graph ΔN vs N for RMD is shown in Fig. 4B. Here we report the values of ΔN sampled at various N in RMD simulations at P = −0.08. The black dashed line, introduced to help the reader to follow the wetting/recovery path, is obtained by connecting the average values of ΔN at each N, 〈ΔN〉N; in the region in which the graph is split into two parts (N > 450) 〈ΔN〉N is computed within the domain of each mode. One notices that the wetting starts (low N) with a symmetric meniscus (ΔN is concentrated about 0). Then, in correspondence of the transition state, the system follows either the bubble-in-the-left (ΔN ≫ 0) or bubble-in-the-right (ΔN ≪ 0) corner paths. With the progress of the wetting process the bubble is absorbed and ΔN tends again to zero. Within the quasi-static approximation, the recovery path is just the reverse of the wetting one, as indicated by the arrows at both ends of the dashed line in Fig. 4B. The ΔN vs. N plots at the different pressures (Fig. SM2†) overlap, indicating that the wetting path is independent of the pressure triggering the process, as predicted by continuum theories of wetting.28
The same analysis is performed on the FFS simulations. In Fig. 5A–C is reported the probability m(ΔN, N) along the wetting at P = −0.005, 0.01 and 0.035. This figure is the FFS analogue of Fig. 4B, representing the probability to observe a value of ΔN at a given level of progress of the wetting as measured by N. One immediately notices that, at a variance with RMD, in this case the wetting path depends on the pressure: the m(ΔN, N) probabilities at different P are sizably different. Like in RMD, at low and moderate liquid pressures (panels A and B) the trajectories are initially (N ≤ 450) symmetric, with m(ΔN, N) concentrated about at ΔN = 0. Then, m(ΔN, N) splits into two branches corresponding to the bubble-in-the-left and bubble-in-the-right corner configurations. At P = 0.035, however, m(ΔN, N) is concentrated about zero all along the process, indicating that at high pressures the mechanism becomes symmetric (Fig. 5C).
Fig. 5 Logarithm of the probability density, log[m(ΔN,N)], along the wetting (A–C) and recovery (D–F) trajectories at different pressures. We refrain from using the same representation of Fig. 4 to highlight that while for RMD one refers to the equilibrium distribution data for FFS obtained by sampling ΔN at prescribed values of N along reactive trajectories. As in the case of Fig. 4, the black dashed line obtained from 〈ΔN〉N is introduced to guide the eye of the reader. At low (P = −0.005) and moderate (P = 0.01) pressures the wetting follows a path consistent with the quasi-static picture of panel D. At higher pressures (P = 0.035) m(ΔN, N) is concentrated about ΔN = 0 all along the wetting path. The recovery always follows a path characterized by initial (N > 400) large positive or negative values of ΔN. However, at small and moderate negative pressures (P = −0.005, −0.08) in the second part of the path (N < 400 and N < 200, respectively) m(ΔN, N) is concentrated about 0, indicating a recovery of the symmetrical morphology of the meniscus. At P = −0.08 we observe multiple transition channels, as highlighted by multiple dashed lines in the central part of the process, between N = 500 and 200. At more negative pressures (P = −0.168) m(ΔN, N) remains centered at large negative or positive values all along the recovery. |
The discrepancies between FFS and RMD suggest that inertia, neglected in the latter approach, plays a crucial role in the wetting of the cavity when the corresponding barrier is low, i.e. when the pressure is high: under these pressure conditions the thermodynamic and viscous forces are insufficient to change the shape of the liquid/gas interface to the minimum free-energy morphology. This claim is supported by the observation that the velocity of advancement of the meniscus, measured by Ṅ, significantly increases with the liquid pressure in correspondence of the transition state. This, indeed, is reflected on the doubling of the Womersley number with pressure, which goes from α2 ∼ 0.035 at P = −0.005 and 0.01 to α2 ∼ 0.07 at P = 0.035.70 In contrast, we observe no significant change of the velocity of bending of the meniscus, estimated by ΔṄ, with P (Fig. 6).
One question arises: to what extent the phenomena we observed on the 10 × 10 (3.5 nm × 3.5 nm) texture are relevant to experimentally accessible lateral scales, ∼10 nm (ref. 15 and 16) or more? As a first, empirical, confirmation that the same phenomena are relevant also on larger scales, we show that the transition from the asymmetric to symmetric path when pressure increases from −0.005 to P = 0.035 is also observed for a 20 × 20 pore (∼7 nm × 7 nm, – see Fig. SI1†). Concerning larger, macroscopic, corrugations, continuum fluid dynamics predicts that the relevance of inertia grows with the characteristic length of corrugation, as shown by the quadratic dependence of the Womersley number α2 = ρL2/μτ on the cavity mouth L. Thus, one expects that the role of inertia in the intrusion mechanism becomes more relevant for larger corrugations.
The present results for the nanoscopic corrugations – the relation between the wetting mechanism, asymmetric or symmetric, and the corresponding barrier – allow us to interpret recent experimental results of Duan and coworkers,24,38 which conflict with well established (macroscopic) quasi-static theories. Present results show that inertia induces a change in the wetting mechanism from asymmetric when the barrier is sizably larger than the thermal energy to symmetric when the barrier is small. At the same time, one expects that in experiments wetting is observed when the barrier is low as under these conditions the wetting time is within the experimental timescale (see eqn (1) and the corresponding discussion). Thus, we conclude that the symmetric wetting mechanism observed in experiments are driven by inertial effects, which are neglected in standard wetting theories. In other words, the inclusion of inertia allows one to reconcile the mismatch between theory and experiments.
These conclusions are qualitatively confirmed considering the transition time τ of FFS wetting simulations (eqn (11)). Indeed, at low and moderate pressures, P = −0.05 and 0.01, corresponding to asymmetric wetting paths, the transition time is 1.5 × 1019 s (∼5 × 1011 years) and 1.5 × 108 s (∼5 years), respectively, well beyond the typical experimental time (microseconds to minutes). In contrast, at high pressure, P = 0.035, at which the transition path is symmetric, the transition time is rather short, 1.5 × 10−8 s (∼15 ns). Transition times in natural units (NU) have been obtained using the standard conversion relation where m is the mass of a water molecule, and ε and σ the Lennard-Jones parameters of the TIP4P2005 water model,71 one of the most widely used water models (see Note72 for additional information). We remark that τNU is purely indicative and one should refrain from performing a quantitative comparison between experimental wetting times and transition times in natural units. There are several reasons for this, among which we recall that while in experiments the system (typically) contains air, our computational sample is completely degassed. This does not affect the wetting path37 but sizably reduces the transition time in simulations with respect to experiments as there is no air that must diffuse into the intruding water for the completion of the process. Here, an important observation is that, consistent with experiments, when the transition time is within the experimental timescale the wetting mechanism is symmetric.
We remark that at the highest pressure, the one at which the transition time is within the experimental timescale, the wetting follows only the symmetric path. This suggests that there is only one wetting channel active at a time, which supports the conclusion of Duan and coworkers24 that there must be an external agent, e.g. accumulation of impurities on the textured walls, inducing the system to follow the asymmetric wetting. Indeed, impurities may pin the meniscus and prevent its symmetric advancing. The effect of impurities accumulating on surface textures on the wetting mechanism will be investigated in detail in a forthcoming study.
We now consider the recovery process by which vapor is formed within the pore and pushes the liquid back to the top of corrugation, restoring the Cassie–Baxter state. We remark that this process is equivalent to cavitation or bubble nucleation under confinement, which are important in many applications.19,73–76 Given the findings on the wetting path, a question naturally arises: do dynamic effects play an important role in the recovery path as well? To address this question we considered three negative liquid pressures (suction), P = −0.005, P = −0.08, and P = −0.168 corresponding to large, intermediate, and negligible free-energy barriers, respectively (Fig. 3 and Table 1). Fig. 5D–F show m(ΔN, N) along reactive recovery trajectories. One notices that in all cases the process is asymmetric and begins with the formation of a bubble in a corner. However, at moderately negative pressures, P = −0.005 and −0.08 – i.e., when the recovery barrier is large – after the transition state the system recovers the symmetric configuration, consistent with the quasi-static picture. Indeed, at P = −0.08 one observes multiple transition channels, as highlighted by multiple dashed lines in the central part of the process, between N = 500 and 200. These multiple transition channels represent different sizes (N) at which the bubble undergoes a transition from the bubble-in-the-corner to a flat morphology. At more negative pressures, instead, the interface becomes even more asymmetric, maintaining this morphology for most of the recovery path. Only when the meniscus pins to one of the corners at the top of the cavity, is the symmetric Cassie–Baxter state recovered. Similar to the wetting case, when the recovery barrier is negligible the velocity of the extruding liquid is too fast for the meniscus to reach the quasi-static configuration at the current level of recovery (Fig. 6). In other words, like in the case of wetting, when the recovery barrier is very low inertial effects dominate. This is confirmed by a sizable growth of the Womersley number at the maximum Ṅ, which passes from ∼0.2 to ∼0.33 for P going from −0.005 to −0.168.
The inertial effects discussed above cause deviations from other equilibrium properties. For example, along wetting and recovery the value of the contact angle and the curvature of the meniscus are different from the ones predicted by the Young and (generalized) Laplace28 equations respectively, describing the shape of a quasi-static liquid/gas interface. For the recovery, these effects are shown in Fig. 7. In FFS simulations at moderately negative pressures (P = −0.005) the contact angle θ at the two points of the meniscus are equal and very close to the Young value, θY = 102°, all along the process. The departure from θY for small vapor bubbles is due to the limited accuracy in the determination of θ under these conditions. In contrast, at P = −0.168, where inertial effects are large, for most of the recovery the values of θ at the two contact points are significantly different, and very different from θY. In particular, we observed a jump in the value of θ at the right contact point when the bubble detaches from the bottom wall (N ∼ 300). Also notice that the curvature of the meniscus when the recovery is almost complete (N ∼ 100) is opposite to the one predicted by the (generalized) Laplace equation.
Contrary to previous quasi-static predictions, it was found that the wetting and recovery transition paths sensitively depend on the liquid pressure, underscoring the importance of liquid inertia. These dynamic effects play a crucial role also in enhancing the difference between the wetting and recovery processes, which are shown to follow symmetric and asymmetric paths, respectively, under typical experimental conditions.
We remark that the revealed phenomenology applies to a variety of physical phenomena well beyond wetting: condensation, cavitation, dynamics of the triple line, micelle formation and many more.
Our new predictions will hopefully stimulate the research community to perform accurate wetting and recovery experiments to identify the mechanism of the processes and their dependence on the thermodynamic conditions and on the direction of the path – wetting or recovery. For example, one could consider working with degassed water, which would help the comparison with simulations, which are typically performed with a pure liquid. Another important aspect to investigate experimentally is the difference between wetting and recovery: to the best of our knowledge, experiments addressing the recovery path are still lacking.
An interesting direction to pursue in the future is to assess the effect of heterogeneity and impurities on the wetting and recovery mechanism. Concerning theory, one could perform simulations on porous systems mimicking the presence of heterogeneities and/or the deposition of impurities on the texture walls.77 Complementarily, one could perform confocal microscopy experiments on systems containing an increasing content of impurities and determine whether the frequency of observation of the asymmetric mechanism has any relation with this.
The implication of present results for simulations is that, when the free-energy barriers characterizing a transition are low, it is crucial to include the dynamics beyond the quasi-static assumption common to many theories (e.g., classical nucleation theory) and rare event methods (e.g., umbrella sampling, restrained MD, string method, and nudged elastic band). From the practical point of view, this will help in developing more accurate strategies to design materials with tailored wetting and recovery properties78 beyond quasi-static theories and methods presently used.
(4) |
(5) |
(6) |
In this work, we have set λ = 0.2, which has already been tested in previous studies to be a good trade-off between the convergence of dΩλ(N*)/dN with λ and the statistical error of the mean force (see, e.g., ref. 21, 55 and 82).
Assuming that the ensemble is at constant temperature, m(r)gλ(N(r) − N*) = exp(−β(V(r) + λ/2(N(r) − N*)2)), which can be sampled by a constant temperature MD driven by the augmented potential Ṽ(r;N*) = V(r) + λ/2(N(r) − N*)2, the so-called Restrained MD. Indeed, one can extend this approach to the isothermal-isobaric ensemble, provided that one uses molecular dynamics suitable to sample this ensemble.
In practice, dΩ(N*)/dN is computed as the time average of λ(N(r) − N*) along the RMD. Indeed, the conditional average of any observable can be computed as the time average of a suitable estimator along an RMD trajectory. We remark that RMD molecular dynamics is only used to sample the conditional ensemble; a single RMD trajectory does not represent the transition and its (fixed) duration bears no relation with the transition time.
(7) |
(8) |
Here, like in previous studies,21,27,28,37,48,52,55,82,83 we use a different approach. We divide the 300000 configurations used to estimate dΩ(N)/dN into M smaller sets from which we obtain the corresponding estimates of the mean force {dΩi(N)/dN}i=1, M. From these, by numerical integration, one obtains a set of free energy curves Ωi(N)1=1, M that can be used to directly compute the variance δΩi(N)2 at each value of N:
(9) |
The FFS algorithm consists of two steps. The first one evolves one (or more) trajectory initialized in the basin A and crossing the first interface N0. From this trajectory one can compute the flux Φ0 of MD trajectories crossing the boundary of the “reactant”. In the second step one evolves many trajectories starting from N0 and reaching either N1 or returning back to N0. From these trajectories one computes the conditional probabilities m(N1|N0) of reaching the next interface. The second step can be iterated for the other interfaces. The total rate of the A → B transition is finally evaluated according to the formula
(10) |
The transition time of the process is the inverse of the rate
τ = 1/kAB. | (11) |
In addition, one obtains a sample of the ensemble of reactive trajectories under the prescribed thermodynamic conditions for the calculation of any statistical quantity along the reactive path.
The first step of the method requires the evolution of a (long) MD trajectory starting from the Cassie–Baxter state. The flux is computed as the number of positive crossings of the first interface, n0, divided by the total duration of the trajectory t
(12) |
The hitting points of this trajectory with the first interface are used to initialize the second step of the algorithm. The hitting points with the second interface are used to initialize further trajectories aiming at the third interface and so forth. The ratio between the successful trajectories, nj+1, those hitting the next interface without returning to N0, and the total number of trajectories fired from the present one, n0j, gives the conditional probability m(Ni+1|Ni) = nj+1/n0j.
To have a statistical consistency, we run a fixed number of trajectories per interface, n0j = 544. Since the number of successful trajectories is lower than the number of fired ones, typically one has an insufficient number of different hitting points from which to start trajectories from the next interface. This problem can be addressed in several ways; here we use an approach in which one starts more than one trajectory from the same configuration. The stochastic dynamics prevents that two trajectories started from the same point overlap.
At a variance with RMD, FFS trajectories represent realistic reactive paths. The duration of these reactive trajectories decreases with the reduction of the barrier (Fig. SI3†). However, this is not the major effect of the barrier; what determines the reduction of the transition time is the increase of the number of reactive trajectories over the total number of attempted jumps from the initial to the final state. In fact, while the duration of the reactive branch of wetting and recovery trajectories, i.e. the part of the trajectory between depinning and complete wetting or vice versa, approximately halves over the range of pressures considered (Fig. SI3†), the transition time, which is the time one has to wait to observe a reactive event, including the time the system spends in the original basin, changes of 27 orders of magnitude.
(13) |
The parameter cij is set to 1 for fluid–fluid and solid–solid particle interactions. For fluid–solid interactions cij is set so that the contact angle θ takes a value close to 100°, typical of silanized or fluorinated surfaces.64–66 The value for cij is chosen following an iterative procedure consisting of three steps:27,55 (i) perform the simulation of a (cylindrical) droplet deposited on a flat surface with a guess value for cij, (ii) determine the value of the contact angle at the present value of cij, and (iii) adjust cij so as to increase/diminish the hydrophobicity of the surface.
The sample for the determination of the contact angle consisted of 54694 fluid and 70000 solid particles. Simulations were run in the NVT ensemble at the same temperature of the RMD and FFS simulations, T = 0.8. The drop (Fig. 8A) has a radius of ∼20. From the MD trajectory one determines the density field (Fig. 8B) and the (Gibbs) liquid/vapor dividing surface. This surface is fitted with a circumference and the contact angle is the derivative of the circumference at its contact point with the nominal position of the surface.
RMD and FFS simulations were run at a constant number of particles, temperature and pressure (NPT ensemble) using the LAMMPS molecular dynamics code.89 Temperature was controlled using the Nosé-Hoover chain thermostat90 with a characteristic time of 0.1. Concerning the pressure, we have adopted the mechanical barostat recently introduced by Marchio et al.,91 which is suited for multi-phase systems. In practice, a solid slab of particles is added above the liquid and an extra force Fext = PA (P is the target pressure and A is the area of the solid slab) is applied to these particles. The time step for the numerical integration of the equation of motion was 0.005.
For RMD simulations we used 35 target values of the order parameter between N* = 160 and N* = 700. At each N* value we performed a 3 × 105 step long RMD simulation.
The number of FFS interfaces is set such that the probability of reaching the next interface is 0.1 < m(Nj+1|Nj) < 0.5. The number of interfaces goes from 15 for the recovery process at very negative pressures to 50 for wetting at moderate pressure. The distance between interfaces varies along the path, going from a difference of 5 particles in the pore close to the transition state to 60 near the product state.
Footnote |
† Electronic supplementary information (ESI) available: Bubble morphology for larger pores and at various pressures. See DOI: 10.1039/c9nr05105h |
This journal is © The Royal Society of Chemistry 2019 |