Activated Wetting of Nanostructured Surfaces: Reaction Coordinates

Subscriber access provided by READING UNIV

Article

Activated Wetting of Nanostructured Surfaces: Reaction Coordinates, Finite Size Effects, and Simulation Pitfalls Matteo Amabili, Simone Meloni, Alberto Giacomello, and Carlo Massimo Casciola J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b07429 • Publication Date (Web): 04 Dec 2017 Downloaded from http://pubs.acs.org on December 4, 2017

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

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

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

The Journal of Physical Chemistry

Activated Wetting of Nanostructured Surfaces: Reaction Coordinates, Finite Size Eects, and Simulation Pitfalls ∗

M. Amabili,

S. Meloni, A. Giacomello, and C.M. Casciola

Dipartimento di Ingegneria Meccanica e Aerospaziale, Università di Roma La Sapienza, Rome, Italy

E-mail: [email protected]

1

ACS Paragon Plus Environment

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

Page 2 of 42

Abstract A liquid in contact with a textured surface can be found in two states, Wenzel and Cassie. In the Wenzel state the liquid completely wets the corrugations while in the Cassie state the liquid is suspended over the corrugations with air or vapor trapped below. The superhydrophobic properties of the Cassie state are exploited for self-cleaning, drag reduction, drug delivery etc., while in the Wenzel state most of these properties are lost; it is therefore of great fundamental and technological interest to investigate the kinetics and mechanism of the Cassie-Wenzel transition.

Computationally, the

Cassie-Wenzel transition is often investigated using enhanced sampling (rare events) techniques based on the use of collective variables (CVs). The choice of the CVs is a crucial task because it aects the free-energy prole, the estimation of the free-energy barriers, and the evaluation of the mechanism of the process.

Here we investigate

possible simulation artifacts introduced by common CVs adopted for the study of the Cassie-Wenzel transition: the average particle density in the corrugation of a textured surface and the coarse-grained density eld at various levels of coarse graining.

We

also investigate possible additional artifacts associated to nite size eects. We focus on a

pillared

surface, a system often used in technological applications. We show that

the use of a highly-coarse-grained density (extreme

coarse graining )

of the uid in the

interpillar region brings to severe artifacts: errors of hundreds of

kB T

in the dierence

of free energy between the Cassie and Wenzel state, of tens of

kB T

in the estimate

of the free-energy barriers, and erroneous wetting mechanisms. A proper description of the wetting mechanism and its energetics apparently requires a ne discretization of the density eld. Concerning the nite-size eects, we have found that the typical systems employed in simulations of the Cassie-Wenzel transition, containing a single pillar within periodic boundary conditions, prevents the complete break of translational symmetry of the liquid-vapor meniscus during the process.

Capturing this break of

symmetry is crucial for describing the transition state along the wetting process, and the early stage of the opposite process the Wenzel-Cassie transition.

2


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


1

Introduction

A liquid in contact with a submerged nanotextured surface can be either in the Cassie 1 or Wenzel 2 state. In the Cassie state, capillary forces maintain a vapor or gas layer below the liquid, which is, thus, in contact only with the top of the surface textures. The presence of a composite vapor-liquid-solid interface makes the surface superhydrophobic, e.g., it reduces the drag 3 or makes the surface self-cleaning. 4 In the Wenzel state, instead, the liquid completely lls the nanotextures with the ensuing loss of all the superhydrophobic properties. In this work we compare dierent atomistic

rare event

simulation approaches to the problem

of the Cassie-Wenzel transition, in order to identify nite-size eects and potential artifacts coming from a coarse-grained description of the process. In normal conditions the transition from Cassie to Wenzel, or the reverse dewetting transition, is a thermally

activated

(stochastic) process, i.e., the system must overcome a free

energy barrier in order to go from one state to the other passing through the transition state. According to the transition state theory, 5 the height of the free energy barrier determines the kinetics of the process. The (average) time needed to observe a transition, or equivalently the lifetime of a given state, depends exponentially on the ratio between the values of the barrier and of the thermal energy kB T , with kB the Boltzmann constant and T the temperature. Transitions times of up to 1200 hours (50 days) have been observed experimentally, 6 which are related to Cassie-Wenzel barriers much larger than kB T . These processes are usually denominated rare events and require specialized techniques to be simulated within reasonable computational times. The Cassie-Wenzel process has frequently been studied by mesoscopic 79 and atomistic 10,11 approaches focusing on the coarse-grained density eld of the uid. One of the objectives of this work is to understand what level of coarse graining is necessary to avoid artifacts signicantly aecting the free energy of the initial, transition, and nal states and the wetting mechanism. The other aspect we want to address is the eect of the size of computational samples on the same aspects, free energetics and mechanism. Indeed, calcu3



lations have been often performed with minimalistic samples including only one unit of the textures (unit cell), 8,1215 which might be insucient to capture qualitative and quantitative aspects of the Cassie-Wenzel transition. Clarifying these aspects is crucial to properly use computational tools to investigate wetting and dewetting of textured surface and to progress towards the

in silico

design of superhydrophobic surfaces with tailored properties.

In this work we use molecular dynamics (MD) in combination with a rare event technique, 16 namely the string method in collective variables, 17 to compute the most probable Cassie-Wenzel transition path and free energetics of a liquid on a surface with nanoscale pillars, a texture type which is often used in technological applications 18 and in experiments; 1921 the nanopillars are also representative of a broad class of surface textures with three-dimensional interconnected vapor domains. The choice of rare-event techniques based on molecular dynamics (MD) is motivated by their fundamental importance: thanks to the minimal assumptions on the multiphase system (solid, liquid, vapor), MD is able to naturally treat phase changes and changes of the topology of interfaces; in addition, it capture eects typical of nanoscale systems, such as density uctuations, layering of the liquid at the walls, and the dependence of the surface tension on the curvature of the meniscus, which are not captured in macroscopic mean-eld models. In other words, the use of MD ensures that the present calculations capture all the essential ingredients of the physics of the system under investigation. The Cassie-Wenzel transition has been studied using several rare event techniques: boxed molecular dynamics, 22 forward ux sampling, 23,24 restrained molecular dynamics, 25,26 (indirect) umbrella sampling, 14 and the string method. 10 All these methods require the use of one or more collective variables (CVs), i.e., a set of observables which can describe the macroscopic conguration of the system along the transition. Because of the computational cost of using several CVs, in most simulations a single CV has been used, namely the number of uid particles within a control volume enclosing the surface corrugations, nl . 14,25,26 Since in the Cassie state, the liquid sits at the top of the textures while in the Wenzel state it wets them 4


Page 4 of 42

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


completely, nl indeed seems a rather natural variable to characterize the process. However, for the specic case of the 2D square cavity it has been shown that a ner coarse-grained density eld is necessary to properly describe the

entire

wetting transition and accurately

determine the corresponding barrier. 10 Consistently, Savoy and Escobedo 23 concluded that the average density of liquid in the corrugations is insucient to describe the Cassie-Wenzel transition of a droplet on a pillar-like surface when the droplet is larger than the inter-pillar distance. Despite the conclusions of these works and the known risk that an unfortunate choice of CVs can critically aect simulation results, 27 a limited attention has been devoted to identify possible artifacts introduced by the use of nl . The conclusions we draw here are not limited by the specic method we use the string in collective variables but have the broader scope of investigating what are suitable collective variables and sample sizes to describe wetting of structured surfaces. We thus expect that analogous conclusions apply to many other rare event techniques in collective variables. 2830 In addition to wetting, the present results might also be relevant to those research themes in which the density eld is used to describe thermally activated events: nucleation of crystals, 3133 protein folding driven by liquid uctuations, 34 membrane poration, 35,36 etc. In these elds there is an interest in developing coarse-grained models, which reduce the computational cost while capturing the relevant physics. The article is organized as follows. In Sec. 2 the atomistic model and the string method are introduced. In Sec. 3 the results are presented and discussed. Finally, Sec. 4 is left for conclusions.

2

Theory

2.1 Atomistic Model The atomistic model used to investigate the Cassie-Wenzel transition consists of a LennardJones (LJ) uid in contact with a nanopillared LJ solid surface (see Fig. 1, left panel). The 5



F

T

Page 6 of 42

F

T

T

F

F

T

T

F

Figure 1: (Left) Atomistic systems formed by 3 × 3 and 1 × 1 pillars. (Right) Discretization of the space corresponding to the dierent sets of CVs for the 3 × 3 (1, 9, and 864 CVs) and 1 × 1 (1 and 96 CVs) pillars systems. The coarse-grained density eld corresponds to the collection of local densities computed within the yellow parallelepipeds shown in the gure which have size δx × δy × δz ≈ 56 × 56 × 29σ , ≈ 18.7 × 18.7 × 29σ , and ≈ 4.6 × 4.6 × 4.8σ for the 1, 9, and 864 CV cases, respectively for the 3 × 3 pillars system. We remark that the 9 and 864 CV cases for the 3 × 3 pillars correspond, in terms of discretization of the space (δx × δy × δz ), to the 1 and 96 CVs for the 1 × 1, respectively. Comparison of results between 1 × 1 and 3 × 3 systems must be restricted to CVs corresponding to equivalent discretizations of the density eld. T and F stand for top and front views, respectively.

6


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


uid is also in contact with a second smooth solid surface of LJ particles at the top of the simulation box, which is used to control the liquid pressure Pl (see below for details). Fluid and solid particles interact via a modied LJ potential

V (rij ) = 4

σ rij

12

σ −c rij

6 !

(1)

where σ and are the characteristic length and energy of the potential, respectively. rij is the interatomic distance between particles i and j . c is the parameter controlling the hydrophobicity 37 of the solid material; here we use c = 0.6, to which corresponds a Young contact angle 38 θY ≈ 113◦ . 25,39 Lennard-Jones units are used throughout the manuscript, i.e., lengths are expressed in

σ and energies in [when it was more convenient to express the (free)-energies in thermal units, kB T , it is explicitly noted]. Conversion in

real units

is obtained using proper values

for and σ , for example /kB = 119.8 K and σ = 3.405 Å for Argon. The uid particles are kept at constant temperature T = 0.8 via a Nosé-Hoover chain thermostat. 40 The pressure of the liquid Pl = 0 is controlled through the application of a constant force Fz on the particles of the upper wall, which acts as a piston, while the bottom wall is kept xed during the simulation. 26 Here we do not use conventional barostats which might introduce additional nite-size artifacts when applied to multiphase systems with domains at dierent pressures, e.g., bubbles, with changing size. 41 Within the quasi-static simulation approach used in the present work, the vapor pressure

Pv depends, essentially, only on T . 42 Thus, the choice of Fz and T determine the driving force of the wetting process, ∆P = Pl − Pv , which is ∆P ≈ 0 in the present case. The two investigated systems are shown in Fig. 1. One consists of a 3 × 3 grid of pillars, and the other of a single pillar, denoted as 1×1 in the following. Periodic boundary conditions are applied in the x and y directions. The pillars are ≈ 18 σ tall and 8 σ thick, and are at a distance of ≈ 11 σ from each other.

7



Page 8 of 42

2.2 Collective variables and the Minimum Free Energy Path Physical processes and chemical reactions in complex environments are conveniently described in terms of collective variables, i.e., a set of nCV observables φ ≡ {φ1 (r), ..., φnCV (r)} that are function of the positions r = {r1 , ..., rnf } of the uid particles. From a macroscopic perspective, the Cassie-Wenzel transition can be described by the advancement of the liquid/vapor meniscus Σlv (sharp interface models) and its twisting and bending. 20,21 An alternative macroscopic approach, the (classical) density functional theory (DFT), 43,44 considers the uid density eld ρ(x) as the fundamental quantity to describe wetting. 7,8,45 It is therefore rather natural to use the coarse-grained density eld as CV in atomistic rare event calculations. In practice, the region of the simulation box comprising the pillars is divided into nCV parallelepipedic cells, and for each cell one counts the number of uid particles inside it: nf 1 X χk (ri ), φk (r) = ∆V i=1

k = 1, ..., nCV ,

(2)

where i runs over the nf uid particles, χk 46 is the (smoothed) characteristic function of the

k -th cell, and ∆V is the volume of the cell. In a quasi-static representation of the wetting process a key quantity is the probability density pφ (N ) that the collective variables φ take on the values N ≡ {N1 , ..., NnCV }. N corresponds to a macroscopic state of the system, i.e., the Cassie, Wenzel or any intermediate state visited during the transition. pφ (N ) is frequently expressed in the logarithmic form and in thermal energy units, the Landau free-energy Ω(N ):

(3)

Ω(N ) = β −1 log pφ (N ) = β −1 log

Z

dr m(r)

nY CV k=1

δ(φk (r) − Nk )

where β = 1/(kB T ) is the inverse thermal energy and δ(·) is the Dirac delta function. The 8


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


probability density function pφ (N ) is expressed in terms of the microscopic measure m(r). High-probability states, i.e. Cassie and Wenzel, are global or local minima of Ω(N ).

pφ (N ) allows to gather information on the mechanism and kinetics of the process. 47,48 In principle, one could run a long MD and compute pφ (N ) from the frequency of the associated events. However, this is not possible when the stable and metastable states of the system are separated by free-energy barriers considerably higher than the thermal energy. Indeed, the time required for the system to visit the relevant regions of the CV space, including the low-probability region corresponding to the transition state, is too large for any present and foreseeable supercomputer. In these cases, one usually resorts to advanced simulation techniques that enhance the sampling of regions of low probability which are relevant for the transition. Here we use the string method in CVs 17 which has already been used in combination with the coarse-grained density eld CV in soft matter systems. 10,34,49 At variance with other methods, the string method does not require to reconstruct the free-energy landscape in the entire CV space. Rather, the objective of the string method is to identify the most probable path for the activated process the wetting path in the present case which is related to a sequence of

snapshots

of density elds of the uid along the Cassie-Wenzel

transition. The uid can follow many wetting paths among which the one identied by the string method is the most probable one. In the presence of large barriers other paths significantly dierent from the one identied by the string method have a negligible probability to be followed by the system. An illustration of multiple possible transition paths for a simple 2D free energy landscape is shown in Fig. 2. The advantage of the string method as compared to sampling techniques aimed at reconstructing the entire free energy landscape is that its computational cost scales only linearly with the number of CVs, while most of the other methods have an exponential dependence on the number of CVs. 50 The string method is, therefore, well suited to deal with the large set of CVs employed in this work. The most probable path {N (α)}, with the path parametrization α ∈ [0, 1], satises the condition 9



ˆ (N (α))∇N Ω(N (α)) M

h

i ⊥

Page 10 of 42

=0

(4)

ˆ (N (α))∇N Ω(N (α)) orthogonal where ⊥ indicates that Eq. 4 refers to the component of M ˆ (N ) is the metric matrix of to the path. ∇N Ω(N ) is the gradient of the free energy and M elements

ˆ ij (N ) = M

Z

dr∇r φi (r) · ∇r φj (r) m(r)

nY CV k=1

δ(φk (r) − Nk ).

(5)

ˆ (N ) is associated to the change of variables from r to φ. 17 α is the independent variable M of the parametric representation of the wetting path N (α). Here we adopt the normalized arc length parametrization of the curve, i.e., α is the fractional length of the arc of the curve between the initial and present points. Thus, in terms of the parameter α, N (0) and

N (1) are the coarse grained densities at the Cassie and Wenzel states, respectively. The physical meaning of Eq. 4 is that along the most probable path N (α) the

eective force

ˆ (N )∇N Ω(N ) has zero components orthogonal to the path. −M

Figure 2: Multiple transition paths connecting the stable and metastable states (absolute and local minima) in an illustrative 2D free energy landscape. The grey line represents the most probable path, the outcome of a string calculation, and the blue dashed lines represent two alternative paths of lower probability. In practice, the path N (α) is discretized in L = 64 images, N (α1 ), . . . , N (αL ), with a constant distance between neighboring images, qP nCV

δN = P L−1 i=1

k=1 (Nk (αi )

qP

− Nk (αi+1 ))2

nCV k=1 (Nk (αi )

10

− Nk (αi+1 ))2


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


(6)

= 1/(L − 1) .

The most probable path is obtained using an iterative algorithm, the string method 51,52 (see Supporting Information), which, starting from a rst guess, produces a path satisfying the

ˆ (N ), can be discretized version of Eq. 4. The main ingredients of Eq. 4, ∇N Ω(N ) and M computed, for example, using Restrained MD (RMD) (Sec 2.3). We remark that the string method is a local search algorithm which converges to the most probable path closest to the rst guess. Several approaches exist to identify other possible paths. The simplest consists in starting the string from dierent initial conditions.

Single ended

approaches related to

the string have been developed to automatically identify possible transition paths starting from a single attractive basin, e.g., the climbing string method. 53 The focus of this work, however, is not nding all possible paths but identifying,

starting from the same rst guess,

the eects on the transition paths and on the free energy connected to dierent CV choices. The rst-guess path is accordingly obtained for all CVs from a high pressure P = 0.6 σ −3 wetting trajectory. The wetting barrier decreases with the pressure and, at this pressure, a wetting event has been observed over a timescale of ∼ 105 timesteps. At each timestep along this spontaneous wetting trajectory, we compute the coarse-grained density eld for all the considered CVs. A polygonal curve passing through these points is constructed and

64 equispaced values (according to Eq. 6) of the coarse-grained density eld are chosen along this polygonal; this is the string of rst guess. This approach, in which the string method is initialized from a continuous atomistic trajectory, prevents potential problems in the calculation of the distance δN between neighboring images due to the translational symmetry of the system: for example, in the simple denition of the distance used here, two congurations corresponding to a translation of the density eld by one or more period in the

x and/or y direction would result into two dierent values of δN , even though they represent the same physical conguration. The correct choice between these two congurations is the one which is the push-forward in time of the previous one. Since, however, the string is a

11



local optimization algorithm (see Supporting Information), the initialization procedure

Page 12 of 42

via

a continuous dynamics described above leads automatically to the selection of the proper density conguration N among its symmetric equivalents. Once the iterative string algorithm is converged and the most likely path is known, one can obtain the associated energetics by computing the line integral of ∇N Ω(N ) along the path, thus obtaining Ω(α). Usually the free energy along the wetting is reported as a function of the liquid fraction Φ(α) = (n(α)−nC )/(nW −nC ), with n(α), nC , and nW the total number of particles in the textures (the region used to dene the CVs of Fig. 1) at the present, Cassie and Wenzel states, respectively. Thanks to the monotonic relation between the fractional arc-length and liquid fraction (see Fig. S3 in the Supporting Information), it is possible to report the free energy in the more convenient and illustrative parametrization Φ. 54 When one uses a single CV, i.e., when the coarse-grained density is replaced by the average density in the relevant domain (the yellow box of Fig. 1), the path is trivial and consists in the increase (or decrease) of the value of the single CV. With a single CV the free-energy landscape in which the most probable path is sought for is 1D; in this case the string method boils down to the simpler RMD. Generally speaking, the choice of CVs (the level of coarse graining of the density eld in the present case) can aect the qualitative and quantitative representation of the transition path. Well established tests exist to validate the quality of the CV set at hand, e.g., the committor test. 17,27 This test consists in numerically computing the probability that trajectories starting from the transition interface reach the products before the reactants and check that the distribution is peaked around 50 %. The transition interface is locally dened as a plane normal to the path passing through the maximum of the free energy (transition state). However, performing the committor test is very expensive. In the case of wetting, due to the slow evolution of the uid toward the Wenzel or Cassie state from the putative transition state, the committor test has been rarely performed. 10 Here we introduce a simpler test allowing to check whether the CV set satises minimal necessary conditions, namely that 12


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


the path does not present unphysical discontinuities, i.e., jumps of the density during the activated process, which obviously cannot be present in a process resulting from continuum and atomistic dynamics.

2.3 Restrained MD and related simulation techniques In order to perform one iteration of the string algorithm one needs to compute ∇N Ω(N )

ˆ (N ) at the set of points {N (αi )}i=1,...,L . In this section we show how these quantities and M can be computed by RMD. Consider Eq. 3 dening the Landau free energy and replace the Dirac delta functions on the r.h.s. by smooth approximations of Gaussian form gλ (φk (r)−Nk ) =

q

2π/(βλ) exp[−βλ/2 (φk (r)−

Nk )2 ]. Within this approximation each component of the gradient of the free energy reads

∂Ωλ (N ) ∂Ω(N ) ≈ = ∂Ni ∂Ni R Q CV gλ (φk (r) − Nk ) dr λ (φi (r) − Ni ) m(r) nk=1 R = QnCV dr m(r) k=1 gλ (φk (r) − Nk ) =

Z

(7)

dr λ(φi (r) − Ni ) p(r|N ) .

Thus, ∂Ω(N )/∂Ni is approximated by the expectation value of the observable λ(φi (r) − Ni ) over the conditional probability density function p(r|N ). Here, we have set λ = 0.2, which is a trade-o between the convergence of ∂Ωλ (N )/∂Ni with λ and the statistical error of the mean force (see Supporting Information Fig. S2(a)). The suitability of this value of λ has also been tested for convergence by analizing the behaviour of Ω with λ (Fig. S2(b)). If m(r) is the measure of the canonical ensemble, m(r) in the form exp (−β(V (r) +

PnCV

k=1

QnCV

k=1

gλ (φk (r) − Nk ) can be cast

λ/2(φk (r) − Nk )2 )), which suggests that the conditional

probability density can be sampled by a constant temperature MD driven by the P CV potential V˜ (r; N ) = V (r) + nk=1 λ/2(φk (r) − Nk )2 , the so-called

Restrained

augmented

MD. Indeed,

this argument can be extended to the isothermal-isobaric ensemble, provided that one uses 13



a molecular dynamics suitable to sample this ensemble. In practice, ∂Ω(N )/∂Ni is computed as the time average of λ(φi (r)−Ni ) along the RMD.

ˆ (N ) can also be computed as a time average along the same RMD, this time considering M the observable ∇r φi (r) · ∇r φj (r). We remark that for a single CV one does not need the string method to compute the wetting/dewetting path, which is just a segment containing the values of the CV describing the process from the initial to the nal state. Thus, in this case RMD allows to compute directly the free energy via the numerical integration of dΩ(N )/dN . In this sense, sometimes we refer to RMD as a proxy for a single-CV string.

3

Results and discussion

This section contains three parts. Sec. 3.1 and Sec. 3.2 focus on the eect of the choice of the CVs (corresponding to dierent levels of coarse graining of the density eld) on the wetting path and on the energetics. Sec. 3.3 compares the results obtained for a corresponding level of coarse graining in the 3 × 3 and 1 × 1 pillars systems.

3.1 3.1.1

3×3

pillars system

Wetting mechanism

Figure 3 shows the evolution of the liquid meniscus along the collapse of the Cassie state for

nCV = 1 (1 × 1 × 1 mesh with a cell of ≈ 56 × 56 × 29 σ ), nCV = 9 (3 × 3 × 1 mesh with cells of ≈ 18.7 × 18.7 × 29 σ ), and nCV = 864 (12 × 12 × 9 mesh with cells of ≈ 4.6 × 4.6 × 4.8 σ - Fig. 1). A color coding is applied to help visualizing the distance of the meniscus from the bottom wall of the surface. The meniscus is identied using the Gibbs dividing surface, i.e., the locus of points where the uid density coincides with the mean value between the liquid ρl and vapor ρv bulk values, (ρl − ρv )/2 ≈ 0.375 σ −3 . Thus, identifying the meniscus from an ensemble of 14


Page 14 of 42

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


A

B

C

1

9 864 Figure 3: Congurations of the meniscus along the wetting process for 1 (red frame), 9 (blue frame), and 864 (green frame) CVs. The color code (bottom left) helps identifying the distance of the meniscus from the bottom wall: blue when the meniscus is at the top of the pillars, red when it touches the bottom, and green in between. In the gure we also report the triple line formed by the liquid/vapor interface with the bottom wall (black) (see Supporting Information for Publication for the movies of the wetting trajectories, les 3x3pillarCV864.gif, 3x3pillarCV9.gif, 3x3pillarCV1.gif). The labels A, B, and C refer to the three steps of wetting process discussed in the text: A is the depinning of the meniscus from the top of the nanopillars and the intrusion of the liquid, B is the touch and spreading of the liquid on the bottom wall, and C is the complete condensation of the rareed liquid. atomistic congurations requires computing a coarse-grained density eld, ρ(x). For this computationally inexpensive post-processing operation we use a ner level of coarse graining than that used for the CVs, namely a 56 × 56 × 29 points grid (cell ≈ 1 × 1 × 1 σ ). The analysis of the meniscus shapes along the wetting process reveals that the collapse mechanism consists of three steps (A-C) (Fig. 3). In step A, the liquid, initially pinned at the top of the pillars, starts to progressively ll the inter-pillar space, with the meniscus assuming an essentially at conformation parallel to the bottom wall. This step is similar for all the CVs sets, and corresponds to the linear increase of the free energy (see also Sec. 3.1.2). In step B, which includes the transition state, the liquid touches and spreads over the bottom wall. The existence of a maximum of the free energy in this domain can be explained using an intuitive argument: i) while the liquid slides down along the pillars with a at meniscus the liquid-solid surface increases resulting in an increase of the energetically unfavorable solid-liquid contribution; ii) after touching the bottom wall the free energy decreases because 15



Page 16 of 42

the two liquid-vapor and vapor-solid interfaces are replaced by a single liquid-solid one. Step B of the wetting path is signicantly dierent in the three cases (Fig. 3). With 1 CV, the liquid touches the bottom wall by completely lling the space between two rows of pillars. Due to the 2D periodicity, this amounts to having an innite strip of liquid touching the bottom of the textured surface parallelled by an innite strip of vapor, with the corresponding meniscus facing the bottom wall at a distance of ca. 6 σ . The process proceeds with the liquid progressively invading neighboring rows or squares of pillars. For both 9 and 864 CVs, step B of the wetting process starts in a similar way, with the liquid touching the bottom between two pillars at a single point. This initial part of the collapse mechanism, with the meniscus touching the bottom substrate between two pillars, is consistent with previous results obtained with atomistic and continuum models. 7,9,11,22,24 However, after the transition state, wetting evolves in a rather dierent way for the two CVs. In the case of 864 CVs the liquid

percolates

in neighboring inter-pillar regions; with 9 CVs

it wets disjoint regions. A remarkable aspect of step B with 1 and 9 CVs is that both present

discontinuities

in the path, corresponding to abrupt changes of the meniscus conguration (Fig. 3). For example, in the case of 1 CV the lling of the corrugations proceeds with the liquid suddenly wetting entire squares or rows among pillars. Concerning the 9 CVs case, comparing the second and third snapshots of step B of the wetting path one notices that the bottom wall between the right-low pair of pillars is initially wet by the liquid, while it gets dry in the next image (see arrows in Fig. 3). As we discuss more extensively below, the origin of these unphysical discontinuities in the paths is the inability of 1 and 9 CVs to distinguish among dierent meniscus congurations: more than one conguration of the meniscus is possible for a given value of the CVs. Step C of the path is similar for all the cases: a rareed liquid region between pairs of pillars transforms into bulk liquid, which eventually completely wets the bottom wall. Videos of the wetting process for the three cases are available in the Supporting Information. 16


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


An important question to address is whether the paths identied with the three sets of CVs are all physically sound. Let us rst revise the concept of transition path in rare event simulations. The path is described as a parametric (discretized) curve in the space of the CVs. The parametric representation is, indeed, only a convenient alternative to the representation of the path as a function of time. 17 Thus, for systems evolving according to continuous dynamics, one would expect that observables change smoothly between successive images of the discretized path. In particular, in the case of wetting of textured surfaces, which is well described by the (classical) density functional theory, 7,15 a necessary condition for the CVs to be suitable to describe the process is that the density eld computed, virtually, on an innitely resolved grid changes in a continuous way along the path. To verify whether this condition is met, we computed the nondimensional Euclidean distance between (ensemble averaged) density elds ρ(x, Φi ) on a very ne 56 × 56 × 29 mesh (cell ≈ 1 × 1 × 1 σ ) at successive images i along the paths, δρ (Φi ), and check whether its value remains (almost) constant. δρ (Φi ) is dened as: qP

δρ (Φi ) =

l,m,n (ρ(xl,m,n , Φi )

ρl

− ρ(xl,m,n , Φi+1 ))2

,

(8)

where the normalization is performed with respect to the bulk liquid density ρl , and the sum P

l,m,n

runs over the cells indices of the very ne mesh, which is the same for all the CVs.

It is observed (Fig. 4) that in steps A and C, δρ (Φi ) ≈ 0.01 independently of the number of CVs, indicating that there are no major dierences in the continuity of the paths identied by the three CVs in these regimes. On the contrary, in step B the paths obtained with 1 and 9 CVs show discontinuities with a sizeable increase of δρ (Φi ), which reaches ≈ 0.15 and

≈ 0.07 in the rst and second case, respectively. These results suggest that 1 and 9 CVs are insucient to describe the wetting process because, in the most important region of the path the one determining the kinetics of the process they are inadequate for describing the continuous trajectory followed by the liquid wetting the surface textures.

17



0.15

1 9 864

0.1

B

δρ

C A

0.05 0

0

0.25

0.5 Φ

0.75

1

Figure 4: Nondimensional Euclidean distance δρ (Φ) between density elds at successive points along the string ρ(x, Φ) (see Eq. 8). Consistently with Fig. 3, red, blue and green curves refer to 1, 9, and 864 CVs, respectively. The labels A, B, and C are dened in Fig. 3 and discussed in the text.

Ω/(kB T )

(a) 1400 900

C B

400

1 9 864

A −100

0

0.25

0.5 Φ

(b)

0.75

1

B

1400

C

A 2000 1500

1150

dΩ/dΦ

Ω/(kB T )

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

Page 18 of 42

B

1000

A

500 0.6

0.7

0.8

Φ

900 0.6

0.8 Φ

1

Figure 5: a) Free-energy proles for 1, 9, and 864 CVs. Error bars are computed via the procedure described in the Supporting Information. b) Zoom of the free-energy prole for 1 CV in region B and C highlighting the change of slope of this curve associated to the morphological transitions discussed in the text. The mean force relative to the rst change of slope at the transition between the A and B is reported in the inset. Inspection of the mean force for 1 CV shows a double discontinuity in this domain. In the Supporting Information the 1 CV free-energy prole and mean force are reported with a ner discretization of the transition path in the B and C domains, conrming the results discussed in the main text. The labels A, B, and C are dened in Fig. 3 and discussed in the text. 18


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


3.1.2

Energetics of the wetting

The dierence in the wetting paths is reected in the dierence among the free energy proles (Fig. 5(a)). In the case of 1 CV we observe an extended, relatively at domain of the free energy prole in step B. Indeed, careful analysis shows that this region is composed of two linear segments with a dierent slope (Fig. 5). The discontinuity of the rst derivative is observed in correspondence of the transitions from one morphology to another (see inset of Fig. 5(b)), namely i) the transition from the

at

meniscus to the liquid partly wetting the

bottom wall (transition from step A to B); ii) the transition from the liquid wetting one row to two rows of pillars; iii) the passage from the liquid wetting only part of the bottom wall, with a layer of

bulk vapor

separating the liquid meniscus from the bottom, to the

liquid wetting the bottom wall, still containing small vapor bubbles between pairs of pillars (transition from step B to C). These three points coincide with the three sharp peaks of δρ . The prole is more regular in the case of 9 and 864 CVs, with a well dened maximum of the free energy between the Cassie and Wenzel states. The quantitative comparison of the three free-energy proles reveals two important aspects. First, the dierence of free energy between the Cassie and Wenzel states, ∆ΩCW , which determines their relative stability, depends on the CVs. 9 and 864 CVs yield consistently

∆ΩCW = 940 kB T , while in the 1 CV case ∆ΩCW is ca. 400 kB T higher. This discrepancy, in turn, aects the barrier of the Wenzel-Cassie transition, reducing its estimate by a corresponding amount in the case of 1 CV. Second, also the Cassie-Wenzel free-energy barrier

∆Ω†CW changes from one CVs set to another. In particular, ∆Ω†CW is approximatively 40 and 80 kB T lower for the 1 and 9 CVs, respectively, as compared to the 864 CVs case. We remark that the dierences in the free energy proles with the three CV sets is not due to an insucient discretization of the path. To show this, we repeated the 1CVs calculation in B and C with three times the number of images (see Supporting Information Fig. S5(b)). The original and more accurate proles show minor quantitative dierences. To further conrm that the dierence in the free energy proles does not depend on an 19



insucient discretization of the path but on the dierence of the CV sets employed, in the region across one of the changes of slope of Ω for we increased the number of point by one order of magnitude (Fig. S5(a)). Also in this case we observed minor quantitative dierences in both the free energy prole and mean force. In the following we explain the origin of the discrepancies of ∆ΩCW and the corresponding barrier using also a simple 2D potential model.

0

25

−1

β

0

−2 −6

−4

−2

0

2

4

x

−2 −6

2 1 0

−10

−4 −2

0 x

2

4 dΩ/dx

50

2

5 3 1 −1

−14 −18 −6

String Marginal tRMD tRMD T = 0

S/kB

1

(b)

α

Ω/kB T

(a) 2 y

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

Page 20 of 42

−3

−4

−2

x −3

−2

0

−1

0

2

1

2

4

x

Figure 6: a) 2D Polynomial potential v(x, y)/(kB T ) = ((0.6x + 0.8y)2 − 1)2 − 2.0(0.8x − 0.6y)(0.6x+0.8y)+0.01(0.8x−0.6y)4 +1.5(0.6x+0.8y)3 in the two arbitrary variables x and y . The potential is characterized by two partly overlapping attractive basins, i.e., depending on the position in y the force along x can be attractive toward either minimum. The color map helps identifying minima and saddle points of the potential. The dashed line represents the string path with two variables; the green line represents the two branches of the discontinuous zero temperature tRMD path computed by applying a restraint on the x variable and starting from basin α. The arrow indicates the point in which the system jumps from one branch to the other. b) Comparison among energy proles obtained with dierent numbers of variables, methods and temperatures. The dashed line is the energy prole along the twovariables string reported as a function of x. The red curve is the exact projection of the 2D R R potential on the variable x, v(x) = −kB T log dy exp[−βv(x, y)]/ dydx exp[−βv(x, y)] and β = 1. The green curve is the free energy obtained by integrating the mean force calculated from a zero temperature tRMD simulation, and the cyan curve is the nite temperature equivalent at kB T = 1. In the inset at the bottom-right of panel b is reported the mean force d¯ v (x)/dx for the exact (red), zero (green) and nite temperature (cyan) tRMD cases. This inset illustrates that in both tRMD cases the forces at the transition between the two branches of the path are discontinuous, with the discontinuity reducing at nite temperature. In the inset at the top-left of panel b the entropy is reported along the nite temperature tRMD path.

Dependence of

∆ΩCW

on the CVs

In general, a direct comparison of the free energies

of dierent sets of collective variables is not possible. In the Supporting Information, Sec. 9, 20


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


we show how this comparison can be made for the specic sets of collective variables used in this work. Among the others, such theoretical arguments indicate that the values of ∆ΩCW should be close for the three CVs while the present results conict with this prediction. The origin of this mismatch can be illustrated considering a model 2D polynomial potential v(x, y) on which the eect of the choice of CVs can be analyzed in detail (Fig. 6(a)). In panel b) the potential along the 2D string path (dashed black line of Fig. 6(a)) is reported; this is the most probable path joining states α and β . We also consider the 1D counterpart of the potential, which is, in our analogy, the equivalent of the free-energy prole obtained with higher coarse graining; this 1D potential is obtained by projecting v(x, y) along the x axis according to Eq. 3: v(x) = −kB T log p(x), where p(x) =

R

dy exp[−βv(x, y)]/ dydx exp[−βv(x, y)] is the R

marginal probability density function. Figure 6(b) shows that ∆vαβ , the energy dierence between the two minima, is indistinguishable if computed along the string path or using the 1D potential v(x); the dierences

along

the path will be analyzed later on when the barriers

are considered. We further compute ∆vαβ as one would obtain from an RMD simulation (thought RMD - tRMD) performed using x as CV; this is expected to provide a numerical approximation to v(x). In tRMD one imagines to start from the minimum α increasing the value of x step by step. In nite time, for each x, the system explores congurations within few kB T from the local minimum y ∗ (x) belonging to the basin of α (green curve in Fig. 6(a)). When the barrier along y separating the α and β basins is suciently small the system jumps into the other basin and successively oscillates around the other branch of the green curve which belongs to the basin β . In Fig. 6(a) the transition from the initial to the nal branch of the path occurs when the barrier along the y direction is 1 kB T . From the ensemble of ν congurations sampled at each x one estimates the mean gradient ν dv(x) 1X ∂v(x, yi ) = dx ν i=1 ∂x

21


(9)


Page 22 of 42

from which the free-energy prole is reconstructed by numerically integrating with respect to x. As noted, the values of ∆vαβ obtained from v(x, y) and the projected v(x) coincide, while ∆vαβ is not captured correctly by tRMD. We remark that the relative error on ∆vαβ estimated via tRMD is ≈ 20%, of the same order of magnitude of the error on ∆ΩCW obtained with 1 CV. The reason of the error on ∆vαβ is that the mean gradient used to compute the free-energy prole with tRMD, Eq. 9, diers from the actual one near the jump between the two valleys of the potential (see the inset of Fig. 6(b). In particular, dv(x)/dx obtained from tRMD is discontinuous when the trajectory passes from one basin to the other. This analysis of the simple 2D potential supports our claim that the dierence in the

∆ΩCW between 1 and 9 or 864 CVs is due to the simulation protocol. The transition from one basin to another in the 2D potential is equivalent to the changes of meniscus morphologies along the wetting process with 1 CV. Thus, we believe that the dierence in

∆ΩCW between 1 and 9 or 864 CVs can be ascribed to the severe discontinuities in the gradient of the free energy in the former case (inset of Fig. 5(b)), which are induced by the extreme coarse graining of the associated density eld in the case of one single CV. 55 Indeed, this problem could be solved combining RMD with parallel tempering, 5658 but, due to the high computational cost, this approach has been rarely adopted in the context of wetting transition. 25

Dependence of the wetting barrier on the CVs

When comparing the free-energy bar-

riers computed via dierent CVs, it is convenient to consider rst the

zero-temperature

limit

for which intuitive arguments hold (discussed in detail below). The subsequent investigation of

nite-temperature eects

leads to important insights into the reliability of these intuitive

arguments and into the actual barrier ordering with number of CVs. Here, the notion of zero temperature

refers to the CVs degrees of freedom, while the atoms, and thus the free-energy

landscape, are always at the physical temperature T = 0.8 /kB . 1 CV is a subset of 9 CVs

22


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


which, in turn, is subset the 864 CVs. Thus, in the present context, the zero-temperature limit of the 1 CV set means that all the remaining 863 (uncontrolled) degrees of freedom take the value corresponding to the local or global minimum at the present value of the restrained CV. An analogous argument holds for the comparison of the

zero temperature

barriers with

1 and 9 CVs, and 9 vs 864 CVs. Let us consider rst the Cassie-Wenzel transition in the limit of zero temperature, and explain the phenomenology using the model 2D polynomial potential. Within this limit and for the x variable only, the system strictly follows the path consisting of the line of the constrained minima of the energy laying in the α basin until the barrier in the orthogonal space is zero (green line in Fig. 6(a)); then the system follows the line of constrained minima in the β basin. The string in 2 variables departs from the zero-temperature tRMD path in the intermediate region, undergoing the transition from the initial to the nal basin at a dierent point. Using a dierent language, the trajectories with one and two variables cross the separatrix, the surface separating the two attractive basins, at dierent points. In Ref. 17 it is shown that the string path crosses the separatrix at the minimum of the potential on this surface, and this point is the actual transition state. Thus, the zero-temperature tRMD transition state has higher or equal energy than the 2D string transition state (see Fig. 6(b)). The argument developed above is not limited to the case of the 2D polynomial potential but is also valid in the case of the CVs for the wetting process discussed in this work (see Supporting Indormation, Sec. 9). One notices that the barriers with 1 and 9 CVs are lower than that with 864, which conicts with the zero-temperature analysis above, which is at the basis of the

intuitive

idea that the barriers should decrease as the number of CVs increases.

This suggests that nite temperature eects play an important role in the ordering of wetting barriers. At nite temperature there are two additional eects to take into account. The rst concerns the point where the transition occurs with a reduced number of CVs. At nite 23



temperature this is attained where the orthogonal barrier is of the order of several kB T , 59 i.e., when the transition time is of the same order of magnitude of the simulation time for computing the mean force. This implies that the transition occurs before than in the zero temperature case, i.e., closer to the true transition state, which results in a reduction of the barrier. The second eect concerns the entropic contribution to the energetics of the transition, which is associated to the number of states the system can explore in the space orthogonal to the collective variables at a given point along the wetting path. The higher is the number of collective variables the lower is the dimensionality of the conguration space the system can explore and, thus, the lower should be the corresponding congurational entropy. For example, in the case of 864 CVs at each point along the string the uctuations allow the density eld to vary over a spatial scale smaller than the cell in which the system is partitioned, i.e. on a spatial scale of ≈ 4.6 × 4.6 × 4.8 σ . In the case of 9 CVs the density can uctuate on a much larger spatial scale of ≈ 18.7 × 18.7 × 29 σ . This eect can be better illustrated in the case of the 2D polynomial potential. At each point along the string in (x, y) the system can take only one conguration, and the entropy is zero. On the contrary, in the case of nite temperature tRMD the system oscillates in the

y direction around the minimum y ∗ (x) at the current value of x; the entropy at each point along the RMD path is dierent from zero. If v(x, y) in the y direction at the present value of x is sti the entropy is low, if it is shallow the entropy is high. In the present case, the potential in the y direction is sti at the initial and nal states and shallow at the transition state. Thus, the entropic contribution, −T ∆S , tends to reduce the barrier as compared to tRMD at T = 0 (see inset of Fig. 6(b). Overall, nite-temperature eects bring the x-only barrier to a value lower than the corresponding string one (cfr. dashed black and solid cyan lines in Fig. 6(b)). We remark that the ordering of the barriers depends on the magnitude of nite-temperature eects and there are no simple

a priori

performing calculations. 24


arguments to estimate it without

Page 24 of 42

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


A similar scenario is likely to occur for the wetting process, in which 1 and 9 CVs might have a higher entropy than 864 at the transition state. For the case of 9 CVs this hypothesis is supported by the visual inspection of the trajectories. The

uctuations

of the meniscus

close to the Cassie state, in which the liquid/vapor interface is pinned at the pillars top, is small in both 9 and 864 CVs cases. On the contrary, at the transition state the meniscus shows larger uctuations for 9 CVs than for 864. To validate this hypothesis we have computed the wetting entropy by subtracting from the free energy the enthalpic term, given by the average Hamiltonian plus the pressure term (SΦ = −1/T [ΩΦ − (hHiΦ + Pl Vl (Φ) + Pv Vv (Φ))], with hHiΦ denoting the ensemble average of the Hamiltonian at the present point of string/RMD and Px and Vx the pressure and volume, respectively, of the phase x, l iquid or

v apor).

It is seen (Fig. 7) that the

entropy of the three CVs initially has a very similar descending trend, which is due to the reduction of the amount of the highly entropic vapor phase (gray dashed line in the gure). In correspondence of the transition states both the 9 and 864 CV cases present an increase of the entropy, which is, however, more pronounced for the 9 CVs. This, indeed, explains why the 9 CVs set has a barrier 80 kB T lower than the 864 one. The scenario is more complex for the 1 CV case. Here, at variance with the other two cases, in step B the entropy further decreases. We believe that this is due to the sudden elimination of a portion of the liquid-vapor interface, which, with its uctuations, contributed to the entropy of the system. At the boundary between B and C, when the 1 CV system changes from a well dened liquid-vapor biphasic system to one made of a liquid plus an extended and diused liquid-solid interface, the entropy suddenly increases. The nal part of the entropy prole is similar for all cases: it decreases following the absorption of the regions of rareed liquid (Fig. 3).

25



1000 −1000 S/kB

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

Page 26 of 42

B

A

−3000

1 9 864 Ideal gas

−5000 −7000

C

0

0.25

0.5 Φ

0.75

1

Figure 7: Entropy along the wetting path for the same sets. The reference is chosen such that the entropy is zero in the Cassie state. The entropy is computed from the free energy by subtracting the enthalpy, which, in turn, is computed as the sum of the expectation value of the Hamiltonian hHi and the −P V term of the liquid and vapor phases. The entropy prole is noisy due to the corresponding noisy signal of hHi. The dashed line represents the entropy of 3 the ideal √ gas computed via the Sackur-Tetrode formula S = kB N (− log(ρv Λ ) + 5/2), where Λ = h/ 3 m kB T is the thermal wavelength, with h the Plank constant and ρv = 0.045σ −3 density of the gas. N is the number of particles in the vapor phase inside the pillars which is computed as N = V (1 − Φ)ρv where V is the volume inside the pillars. 3.1.3

Summary

Let us close this section drawing some conclusions on the eect of the degree of coarse graining on the study of the wetting process. Extreme coarse-graining leads to discontinuities in the wetting path. This results in severe errors in the estimation of ∆ΩCW , up to 400 kB T , and the wetting barrier, which is underestimated by up to ≈ 80 kB T as compared to the ner density eld case. Overall, this suggests that quantitative predictions require a careful choice of the CVs. We remark that these eects are not due to the special method used RMD or string: analogous artifacts would be observed also with US, BXD, TAMD, or other techniques. A detailed analysis of string/RMD

vs

US and BXD is discussed in Supporting

Information.

3.2

1×1

pillar system

For the 1 × 1 pillar system we perform an analysis similar to the one of Sec. 3.1. In this case we considered 1 and 96 CVs, having the same grid spacing of the 9 and 864 CVs cases of the 3 × 3 pillars, respectively. Fig. 8 reports ρ(x, Φi ), Ω(Φi ), and δρ (Φi ). The general 26


Page 27 of 42

B

A

C

(a) 1

96

(b)

A

0.03 δρ

B

170

(c)

0.04

150

C

0.02

0

0.25

0.5

0.75

Ω/(kB T )

100

1

Φ

B

50

A

0 0

C

150

0.01

Ω/(kB T )

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


0.25

0.5 Φ

1 96 0.75

130

B

110

C

90

1

1 96

0.9

0.95 Φ

1

Figure 8: a) Liquid meniscus along the Cassie-Wenzel transition for the 1 × 1 system. Black and violet frames correspond to 1 and 96 CVs, respectively. The color code is the same as in Fig. 3: blue when the meniscus is at the top of the pillar, red when it touches the bottom, and green in between (see Supporting Information for Publication for the movies of the wetting trajectories, les 1x1pillarCV96.gif, 1x1pillarCV1.gif, and 1x1regionC.gif). b) Free-energy proles for the two sets of CVs. The inset shows δρ as a function of Φ. c) Magnication of the region C of the free-energy prole together with the c(1 − Φ)2 (grey dotted) and a(1−Φ)1/2 +b(1−Φ) (grey dashed) curves discussed in the text. The coecients a, b, and c are obtained by tting the two functions in the domain C. To obtain a more careful comparison of the 1 and 96 CVs proles in this domain, the integration has been performed right to left, which results in lower error bars.

27



characteristics of the wetting mechanism observed for the 3 × 3 pillars system are preserved also in the smaller sample: the liquid initially depins from the pillars top, slides along their side with a at meniscus (A), then touches the bottom wall in a point between a pair of pillars (B), and nally the vapor bubble is absorbed (C). It is worth remarking that an analogous mechanism has been reported in the literature for the case of water wetting a 1 × 1 pillared surface. In this case, simulations were performed using INDUS with 1 CV analogous to the one used in the present work. 14 The 96 CVs case presents a smooth transition from one regime to another, with δρ remaining almost constant at 0.01. On the contrary, with 1 CV

δρ shows abrupt jumps in correspondence of the morphological transitions at the boundaries between domains A and B and B and C. Videos of the wetting process for the two cases are available in the Supporting Information. The 1 and 96 CVs free energy proles present no qualitative dierences: both curves are characterized by a single, well dened maximum in correspondence of the conguration in which the meniscus touches the bottom wall. Concerning quantitative aspects, while ∆ΩCW is the same in both cases, the barrier of 1 CV is 20 kB T lower than for 96 CVs. This is due to the nite temperature eects explained in the previous section. For the 1 × 1 pillar case we also investigated artifacts in the early stage of the dewetting mechanism. To make our analysis more accurate, the path with 96 CVs in the range

Φ ∈ [0.9, 1] has been discretized in 32 images (Fig. 8(c)). In this region, the curvature of the free energy proles with 96 and 1 CVs presents sizable dierences. With 1 CV a parabolic prole extending over a broad range of lling fraction is observed, which is generally attributed to Gaussian uctuations of the liquid density as discussed in the Lum-ChandlerWeeks theory of hydrophobicity. 60 A similar trend has been observed by other authors using similar CVs. 14,25,26,61 On the contrary, with 96 CVs the parabolic trend is observed over a much narrower domain, after which the free energy prole shows a dependence of the type

a(1−Φ)1/2 +b(1−Φ) as predicted by macroscopic theories, which accounts for the liquid-vapor and solid-vapor interface energies 62 (see Supporting Information for further details). This is 28


Page 28 of 42

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


reected on an earlier formation of the liquid-vapor meniscus, as shown in the snapshots of the C region of Fig. 8(a). The dewetting mechanism identied with 1 CVs is characterized by a strong initial (≈ 8 %) depletion of the density of the liquid in the corrugation before the meniscus is formed. On the contrary, the path obtained with 96 CVs indicates a more modest ≈ 3 − 4 % density depletion before the liquid-vapor interface forms. 63 Considering that the 96 CVs is a superset of 1 CV (overall density), i.e., the 96 CV-space includes and extends the 1 CV one, the dierent behavior of the free energy indicates that the overall density is inadequate for describing the dewetting path. In several works 14,26,61 the use of the overall density as the only descriptor of the wetting/dewetting process has been justied on the basis of the Lum-Chandler-Weeks theory 60 of hydrophobicity. However, the above observation that 1 CV is insucient to characterize thermally activated wetting/dewetting processes suggests that the extension of the theory of the hydrophobic eect out of its original equilibrium scope may bring to overlooking important aspects of the transition. Indeed, uctuations of the overall density happen during dewetting but for the process to take place they must be accompanied by other events, e.g., the formation, bending, and displacement of the liquid meniscus.

3.3 Size eect: comparison between 3 × 3 and 1 × 1 pillars systems In the previous sections we focused on the eect of the choice of CVs on the wetting mechanism, energetics and kinetics; here we concentrate on the eect of the size of the surface sample, the number of pillars in the x and y directions. To this end we compare results obtained with 96 CVs for the 1 × 1 surface and 864 CVs for the 3 × 3 one. At the transition state the two systems present apparently similar congurations. In both cases the liquid touches the bottom wall at a point between two pillars (Fig. 9). However, the 1 × 1 surface cannot reveal the complex nature of the nal part of the wetting. Indeed, for the 3 × 3 case we observe the formation of a percolating (random) network of vapor bubbles, which is absorbed when the system approaches the Wenzel state. 11 At variance with this scenario, 29



for the 1 × 1 surface, due to the constraints imposed by the periodic boundary conditions, the network of vapor bubbles forms between a pillar and its periodic image along one of the two

lattice

directions, x and y , 64 which brings to the formation of unphysical liquid stripes

(compare the top and bottom panels of Fig. 9).

Figure 9: a) Meniscus at the transition states of the 1 × 1 (center and top) and 3 × 3 (bottom) systems, for which we considered 96 and 864 CVs, respectively, i.e., they have the same grid spacing. The (black) three-phase contact line at the bottom wall helps comparing the morphology of the menisci at the transition state. This gure clearly illustrates that the 1 × 1 system is insucient for identifying possible breaks of the translational symmetry of the meniscus. b) Menisci at selected values of liquid fraction in the C domain. These snapshots show that the 1 × 1 system is insucient for capturing the percolating network of vapor bubbles characterizing this part of the wetting/dewetting path. Concerning the quantitative comparison of results, Fig. 10 reports the free-energy proles for the 3 × 3 and 1 × 1 pillars systems. Given the dierent size of the two systems, such an analysis is perfomed multiplying by a factor 9 the free-energy prole of Fig.8b. While, as expected, the free-energy dierence between Cassie and Wenzel states, ∆ΩCW , are in agreement, the barrier of the 1 × 1 pillar system is ≈ 200 kB T higher than the 3 × 3 one. Moreover, also the position of the transition states is dierent in the two cases: Φ ≈ 0.75 for 30


Page 30 of 42

Page 31 of 42

the 3 × 3 pillars system as compared to Φ ≈ 0.80 for the 1 × 1 one. This is the result of small size of the sample, more specically, of the 9 (periodic) repetition of the bent meniscus in the replicated 1 × 1 system, see Fig. 9. A fairer comparison can be obtained considering the free energy of the transition state plus eight times that of the last at meniscus conguration of the 1 × 1 system. However, also in this case the 1 × 1 barrier is 180 kB T higher than the 3 × 3. This brings us to the conclusion that the 1 × 1 system is aected by important nite-size errors, due to the enforcement of an articial symmetry. 1400 Ω/(kB T )

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


900

B C

400 3x3 1x1

A −100

0

0.25

0.5 Φ

0.75

1

Figure 10: Free-energy prole for the 1 × 1 (96 CV, violet line) and 3 × 3 (864 CV, green line) systems. The 1 × 1 prole is obtained multiplying by 9 the free-energy prole of Fig. 8(b).

4

Conclusions

In this work the Cassie-Wenzel transition has been studied using the string method in collective variables. The collective variable used is the density eld at various levels of coarse graining. Results show the importance of the choice of the collective variables in the case of wetting and dewetting: extreme coarse graining of the density eld, often used in the recent literature, introduces qualitative and quantitative artifacts. These artifacts include unphysical discontinuities in the Cassie-Wenzel transition path, incorrect estimates of the free-energy dierence between the Cassie and Wenzel states and of the associated barrier. We remark that these artifacts are not due to the special rare-event technique used in the present work: analogous eects would have been observed with any other method based on the sampling of the collective variable space. 31



Page 32 of 42

We have also investigated the eect of the size of the sample representing the surface, namely the number of pillars in the simulation box. Results show that a single unit cell (1×1 pillar system) is insucient when dealing with interconnected textures. Indeed, the periodic boundary conditions impose a non-physical translational symmetry of the meniscus leading to an incomplete picture of the wetting mechanism. At least for pillars, in order to capture the detailed wetting mechanism, the simulated system must be large enough to capture the local deformation of the meniscus forming a

liquid nger

which touches the bottom wall

between two pillars. Concerning dewetting, the branch between the Wenzel and transition states consists of a (random) percolating network of vapor bubbles, which cannot be formed in the 1 × 1 pillar system. These simulation pitfalls preclude the use of too small samples for deriving design criteria for preventing the wetting or improving the superhydrophobicity recovery. The results of this work seem to indicate that the simplest simulations of thermally activated wetting with few collective variables and small periodic boxes should be taken cum grano salis :

the information they yield is only qualitative and could possibly be used as a

computationally inexpensive way to identify trends or to perform parametric studies. The detailed physics of the wetting/dewetting processes could be quantitatively dierent from the simulated one and the value of the free-energy barriers could be o by tens of kB T . Depending on the application and on the goal of the work, one should choose what is the best compromise between a physically sound picture of the phenomena and computationally convenient ways to explore dierent cases in order to devise design strategies.

A posteriori,

one could evaluate the choice of CVs: typical symptoms of the deciency of the adopted description is the presence of sudden jumps of the slope of the free-energy proles or the abrupt change of the meniscus morphology, as numerically revealed by peaks in δρ (Φi ).

32


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


Acknowledgement The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement n. [339446]. We acknowledge PRACE for awarding us access to resource FERMI based in Italy at Casalecchio di Reno.

Supporting Information Available See Supporting Information for additional simulations data, study of the convergence of the string and RMD methods, and detail on the macroscopic estimation used to t the data in Fig. 8c.

References (1) Cassie, A.; Baxter, S. Wettability of porous surfaces.

Trans. Faraday Soc. 1944, 40,

546551. (2) Wenzel, R. N. Resistance of solid surfaces to wetting by water. 28,

Ind. Eng. Chem. 1936,

988994.

(3) Rothstein, J. P. Slip on superhydrophobic surfaces.

Annu. Rev. Fluid Mech. 2010, 42,

89109. (4) Bhushan, B.; Jung, Y. C.; Koch, K. Self-cleaning eciency of articial superhydrophobic surfaces.

Langmuir 2009, 25,

32403248.

(5) Eyring, H. The activated complex in chemical reactions. 115.

33


J. Chem. Phys. 1935, 3,

107


Page 34 of 42

(6) Xu, M.; Sun, G.; Kim, C.-J. Innite lifetime of underwater superhydrophobic states. Phys. Rev. Lett. 2014, 113,

136103.

(7) Ren, W. Wetting transition on patterned surfaces: transition states and energy barriers. Langmuir 2014, 30,

28792885.

(8) Panter, J. R.; Kusumaatmaja, H. The Impact of Surface Geometry, Cavitation, and Condensation on Wetting Transitions: Posts and Reentrant Structures. J. Phys.: Matter 2017, 29,

Cond.

084001.

(9) Pashos, G.; Kokkoris, G.; Papathanasiou, A.; Boudouvis, A. Wetting transitions on patterned surfaces with diuse interaction potentials embedded in a Young-Laplace formulation.

The Journal of chemical physics 2016, 144,

034105.

(10) Giacomello, A.; Meloni, S.; Müller, M.; Casciola, C. M. Mechanism of the Cassie-Wenzel transition via the atomistic and continuum string methods.

J. Chem. Phys. 2015, 142,

104701. (11) Amabili, M.; Giacomello, A.; Meloni, S.; Casciola, C. M. Collapse of superhydrophobicity on nanopillared surfaces.

Phys. Rev. Fluids 2017, 2,

034202.

(12) Hemeda, A.; Tafreshi, H. V. General formulations for predicting longevity of submerged superhydrophobic surfaces composed of pores or posts.

Langmuir 2014, 30,

10317

10327. (13) Zhang, Z.; Kim, H.; Ha, M. Y.; Jang, J. Molecular dynamics study on the wettability of a hydrophobic surface textured with nanoscale pillars. 2014, 16,

Phys. Chem. Chem. Phys.

56135621.

(14) Prakash, S.; Xi, E.; Patel, A. J. Spontaneous recovery of superhydrophobicity on nanotextured surfaces.

P. Natl. Acad. Sci. USA 2016, 113,

34


55085513.

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


(15) Tretyakov, N.; Papadopoulos, P.; Vollmer, D.; Butt, H.-J.; Dünweg, B.; Daoulas, K. C. The Cassie-Wenzel transition of uids on nanostructured substrates: Macroscopic force balance versus microscopic density-functional theory.

J. Chem. Phys. 2016, 145,

134703. (16) Bonella, S.; Meloni, S.; Ciccotti, G. Theory and methods for rare events. B 2012, 85,

Eur. Phys. J.

119.

(17) Maragliano, L.; Fischer, A.; Vanden-Eijnden, E.; Ciccotti, G. String method in collective variables: minimum free energy paths and isocommittor surfaces. 2006, 125,

J. Chem. Phys.

024106.

(18) Verho, T.; Korhonen, J. T.; Sainiemi, L.; Jokinen, V.; Bower, C.; Franze, K.; Franssila, S.; Andrew, P.; Ikkala, O.; Ras, R. H. Reversible switching between superhydrophobic states on a hierarchically structured surface. 2012, 109,

P. Natl. Acad. Sci. USA

1021010213.

(19) Checco, A.; Ocko, B. M.; Rahman, A.; Black, C. T.; Tasinkevych, M.; Giacomello, A.; Dietrich, S. Collapse and reversibility of the superhydrophobic state on nanotextured surfaces.

Phys. Rev. Lett. 2014, 112,

216101.

(20) Papadopoulos, P.; Mammen, L.; Deng, X.; Vollmer, D.; Butt, H.-J. How superhydrophobicity breaks down.

P. Natl. Acad. Sci. USA 2013, 110,

32543258.

(21) Lv, P.; Xue, Y.; Liu, H.; Shi, Y.; Xi, P.; Lin, H.; Duan, H. Symmetric and Asymmetric Meniscus Collapse in Wetting Transition on Submerged Structured Surfaces. Langmuir 2014, 31,

12481254.

(22) Savoy, E. S.; Escobedo, F. A. Simulation study of free-energy barriers in the wetting transition of an oily uid on a rough surface with reentrant geometry. 28,

1608016090.

35


Langmuir 2012,


Page 36 of 42

(23) Savoy, E. S.; Escobedo, F. A. Molecular simulations of wetting of a rough surface by an oily uid: Eect of topology, chemistry, and droplet size on wetting transition rates. Langmuir 2012, 28,

34123419.

(24) Shahraz, A.; Borhan, A.; Fichthorn, K. A. Kinetics of droplet wetting mode transitions on grooved surfaces: Forward ux sampling.

Langmuir 2014, 30,

1544215450.

(25) Giacomello, A.; Meloni, S.; Chinappi, M.; Casciola, C. M. CassieBaxter and Wenzel states on a nanostructured surface: phase diagram, metastabilities, and transition mechanism by atomistic free energy calculations.

Langmuir 2012, 28,

1076410772.

(26) Amabili, M.; Lisi, E.; Giacomello, A.; Casciola, C. M. Wetting and cavitation pathways on nanodecorated surfaces.

Soft matter 2016, 12,

30463055.

(27) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Transition path sampling: Throwing ropes over rough mountain passes, in the dark. 2002, 53,

Annu. Rev. Phys. Chem.

291318.

(28) Torrie, G. M.; Valleau, J. P. Nonphysical sampling distributions in Monte Carlo freeenergy estimation: Umbrella sampling.

J. Comp. Physics. 1977, 23,

187199.

(29) Patel, A. J.; Varilly, P.; Chandler, D. Fluctuations of Water near Extended Hydrophobic and Hydrophilic Surfaces.

J. Phys. Chem. B 2010, 114,

16321637.

(30) Glowacki, D. R.; Paci, E.; Shalashilin, D. V. Boxed Molecular Dynamics: A Simple and General Technique for Accelerating Rare Event Kinetics and Mapping Free Energy in Large Molecular Systems.

J. Phys. Chem. B 2009, 113,

(31) Kelton, K. F.; Greer, A. L. and biology ;

1660316611.

Nucleation in condensed matter: applications in materials

Pergamon, 2010.

(32) Page, A. J.; Sear, R. P. Heterogeneous nucleation in and out of pores. 2006, 97,

065701. 36


Phys. Rev. Lett.

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


(33) Page, A. J.; Sear, R. P. Crystallization controlled by the geometry of a surface. Chem. Soc. 2009, 131,

J. Am.

1755017551.

(34) Miller, T. F.; Vanden-Eijnden, E.; Chandler, D. Solvent coarse-graining and the string method applied to the hydrophobic collapse of a hydrated chain. USA 2007, 104,

P. Natl. Acad. Sci.

1455914564.

(35) Fuhrmans, M.; Marelli, G.; Smirnova, Y. G.; Müller, M. Mechanics of membrane fusion/pore formation.

Chem. phys. lipids 2015, 185,

109128.

(36) Smirnova, Y. G.; Fuhrmans, M.; Vidal, I. A. B.; Müller, M. Free-energy calculation methods for collective phenomena in membranes.

J. Phys. D Appl. Phys. 2015, 48,

343001. (37) Cottin-Bizonne, C.; Barrat, J.-L.; Bocquet, L.; Charlaix, E. Low-friction ows of liquid at nanopatterned interfaces.

Nat. Mater. 2003, 2,

237240.

(38) The Young contact angle is the angle formed between a at surface and the tangent to a sessile liquid droplet deposited on it at the liquid/solid contact point. (39) Amabili, M.; Giacomello, A.; Meloni, S.; Casciola, C. Unraveling the Salvinia Paradox: Design Principles for Submerged Superhydrophobicity.

Adv. Mater. Interf. 2015, 2,

1500248. (40) Martyna, G. J.; Klein, M. L.; Tuckerman, M. NoséHoover chains: the canonical ensemble via continuous dynamics.

J. Chem. Phys. 1992, 97,

26352643.

(41) Marchio, S.; Meloni, S.; Giacomello, A.; Valeriani, C.; Casciola, C. M. Pressure Control in Interfacial Systems: Atomistic Simulations of Vapor Nucleation.

In preparation

(42) Amabili, M.; Giacomello, A.; Meloni, S.; Casciola, C. Intrusion and extrusion of a liquid on nanostructured surfaces.

J. Phys: Condens. Mat. 2016, 29,

37


014003.


(43) Evans, R. In

Fundamentals of Inhomogeneous Fluids ;

Page 38 of 42

Henderson, D., Ed.; Marcel

Dekker: New York, 1992; Chapter 3, pp 85176. (44) Löwen, H. Density functional theory of inhomogeneous classical uids: recent developments and new perspectives.

J. Phys.: Condens. Matter 2002, 14,

11897.

(45) Pashos, G.; Kokkoris, G.; Boudouvis, A. G. Minimum energy paths of wetting transitions on grooved surfaces.

Langmuir 2015, 31,

30593068.

(46) χk (r), the characteristic function, is equal to 1 if r is in k -th cell and 0 otherwise. χk (r) can be expressed as the product of two Heaviside step functions per Cartesian direction:

χk (r) =

k,b k,e k,b and Rlk,e the b egin and e nd l=x,y,z Θ(rl − Rl ) 1 − Θ(rl − Rl ) , with Rl

Q

of the k-th cell in the l direction, and rl the component of the particle's position in the same direction. In practice, the Heaviside step function is replaced by a smooth approximation, in this case, a Fermi function. We have found this choice, which is consistent with the literature, 34 convenient, but others are possible. (47) Indeed, to obtain a complete understanding of the transition mechanism and kinetics one needs also the committor function qφ (N ), i.e. the probability that a trajectory in

φ(r) = N reaches the nal state rst. In the string method discussed in this section, which is based on the backward Kolmogorov equation, the committor function is taken into account implicitly. (48) Vanden-Eijnden, E.

Computer Simulations in Condensed Matter Systems: From Mate-

rials to Chemical Biology Volume 1 ;

Springer, 2006; pp 453493.

(49) Müller, M.; Smirnova, Y. G.; Marelli, G.; Fuhrmans, M.; Shi, A.-C. Transition path from two apposed membranes to a stalk obtained by a combination of particle simulations and string method.


228103.

(50) Laio, A.; Rodriguez-Fortea, A.; Gervasio, F. L.; Ceccarelli, M.; Parrinello, M. Assessing the accuracy of metadynamics.

J. Phys. Chem. B 2005, 109,

38


67146721.

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


(51) Weinan, E.; Ren, W.; Vanden-Eijnden, E. String method for the study of rare events. Phys. Rev. B 2002, 66,

052301.

(52) E, W.; Ren, W.; Vanden-Eijnden, E. Simplied and improved string method for computing the minimum energy paths in barrier-crossing events. 126,

J. Chem. Phys. 2007,

164103.

(53) Ren, W.; Vanden-Eijnden, E. A climbing string method for saddle point search. Journal of chemical physics 2013, 138,

The

134105.

(54) We remark that the passage from the coarse-grained density eld φ to the liquid volume fraction in the corrugation Φ is a non-invertible map, i.e., one has one and only one value of Φ for each value of the density eld φ but the opposite that there is only one value of the density eld φ corresponding to the liquid volume fraction Φ is false. In other words, the fact that there is a linear relation between the fractional arc-length, and then the associated density eld φ, and Φ does not mean that Φ and φ, i.e., CVs at dierent level of coarse graining, are equivalent for describing the wetting process. Indeed, the observation that Φ grows monotonically with the fractional arc-length just means that the wetting process occurs with a continuous increase of volume of the liquid in the cavity, a fact that is not obvious but intuitive. (55) This problem is also at the origin of the hysteresis observed with other rare event techniques when one uses insucient or inadequate CVs. (56) Swendsen, R.; Wang, J. Replica Monte Carlo simulation of spin glasses. Lett. 1986, 57,

Phys. Rev.

26072609.

(57) Orlandini, S.; Meloni, S.; Ciccotti, G. Combining Rare Events Techniques: Phase Change in Si Nanoparticles.

J. Stat. Phys. 2011, 145,

812830.

(58) Orlandini, S.; Meloni, S.; Colombo, L. Order-disorder phase change in embedded Si nanoparticles.

Phys. Rev. B 2011, 83,

235303. 39



Page 40 of 42

(59) The transition from the initial to the nal branch of the path is performed when the barrier along the y direction is zero and 1kB T for T = 0 tRMD and tRMD, respectively. The rst condition is met at the green arrow of Fig. 6(a). The second condition occurs very close to this point. (60) Lum, K.; Chandler, D.; Weeks, J. D. Hydrophobicity at small and large length scales. J. Phys. Chem. B 1999, 103,

45704577.

(61) Remsing, R. C.; Xi, E.; Vembanur, S.; Sharma, S.; Debenedetti, P. G.; Garde, S.; Patel, A. J. Pathways to dewetting in hydrophobic connement. USA 2015, 112,

P. Natl. Acad. Sci.

81818186.

(62) Giacomello, A.; Chinappi, M.; Meloni, S.; Casciola, C. M. Metastable wetting on superhydrophobic surfaces: continuum and atomistic views of the Cassie-BaxterWenzel transition.


226102.

(63) We remark that similar results are present also in the 3 × 3 pillar case, in which a dierent slope of Ω in the domain C is apparent (Fig. 5). The detailed analysis of the early stage of dewetting, which requires a ne discretization of the path with a large number of string images, was performed for the 1 × 1 case which is computationally more convenient. (64) We remark that another equivalent wetting path is possible, i.e., one in which the liquid wets the surface along the perpendicular direction. Our string simulations show only the one with the symmetry selected by the rst guess used to initialize the string.

40


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


Graphical TOC Entry Coarse-graining levels

Finite-size effects

41


Coarse-graining levels Finite-size effects The Journal of Physical Chemistry Page 42 of 42

1 2 3 4 5 6


Activated Wetting of Nanostructured Surfaces: Reaction Coordinates

Recommend Documents