Activated wetting of nanostructured surfaces: reaction coordinates, finite size effects, and simulation pitfalls
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 affects the free-energy profile, 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 field at various levels of coarse graining. We also investigate possible additional artifacts associated to finite size effects. We focus on a pillared surface, a system often used in technological applications. We show that the use of the average density of fluid in the interpillar region brings to severe artifacts in the relative Cassie-Wenzel stability and in the transition barrier, with an error on the estimates of the free energies of up to hundreds of and wrong wetting mechanisms. A proper description of the wetting mechanism and its energetics requires a rather fine discretization of the density field. Concerning the finite size effects, we have found that the typical system employed in the Cassie-Wenzel transition, containing a single pillar within periodic boundary conditions, prevents the 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.
A liquid in contact with a submerged nanotextured surface can be either in the Cassie cassie1944 or Wenzel wenzel1936 state. In the Cassie state, a vapor or gas layer sustains 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, i.e. characterized by properties such as enhanced slippage rothstein2010 or self-cleaning. bhushan2009 In the Wenzel state, instead, the liquid completely fills the nanotextures with the ensuing loss of all the superhydrophobic properties.
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. According to the transition state theory eyring1935, the height of the free energy barrier determines the kinetics of the process. Indeed, 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 , with the Boltzmann constant and the temperature. The presence of large free-energy barriers has been confirmed experimentally. Indeed, very long Cassie state lifetimes, of the order of tens of minutes to tens of hours, have been observed for textured surfaces subject to moderate liquid pressures ( higher than the atmospheric pressure xu2014; bobji2009).
Molecular dynamics, thanks to its minimal assumptions on the model of the three-phase system (solid, liquid, vapor), seems an optimal tool to investigate the mechanism and kinetics of the Cassie-Wenzel transition on nanotextured systems, and their dependence on the surface topography and on the thermodynamic conditions. However, due to the long transition time connected with large free energy barriers, special simulation techniques are needed to investigate the collapse of the superhydrophobic state. Thus, the Cassie-Wenzel transition has been studied using several advanced sampling techniques: boxed molecular dynamics, savoy2012b forward flux sampling, savoy2012a restrained molecular dynamics, giacomello2012; amabili2016 (indirect) umbrella sampling, prakash2016 and the string method. giacomello2015 All these methods require the use of one or more (typically few) collective variables (CVs), i.e., a set of observables which can describe the macroscopic configuration of the system along the transition. In most of the simulations a single CV has been used, namely the number of fluid particles within a control volume enclosing the surface corrugations, . giacomello2012; savoy2012a; amabili2016; prakash2016 In a two-phase liquid-vapor system this CV is related to the amount of liquid in the same region. Since in the Cassie state the liquid is pinned at the top of the texture while in the Wenzel state it wets the corrugations completely, seems the most natural variable to characterize the process. has several advantages: i) it is relatively simple to use in atomistic simulations; ii) it is an atomistic equivalent of the volume of liquid in the corrugations, which allows to establish a connection and compare results with the continuum theory of wetting of confined systems, CREaM giacomello2012a; iii) considering the reverse process, dewetting, it is related to the liquid density and allows to establish connections with the Lum-Chandler-Weeks’s theory of
hydrophobicity lum1999 (see Refs. remsing2015; amabili2016). However, for the specific case of the 2D square cavity it has been shown that the (coarse-grained) density field is necessary to properly describe the wetting transition and accurately determine the corresponding barrier. giacomello2015 Despite the conclusions of this latter work, and the known risk that an unfortunate choice of CVs can critically affect simulation results, bolhuis2002 a limited attention has been devoted to identify possible artifacts introduced by the use of . This is especially important in the case of surfaces with 3D interconnected cavities, such as pillared structures, for which it has been reported that the transition takes place with changes in the morphology of the liquid-vapor meniscus. prakash2016
Another aspect that might introduce artifacts in the study of the wetting transition is the size of the simulation box. In the case of artificial textured surfaces, this is translated into the number of unit cells which are included in the computational sample, e.g., the number of pillars for a pillared surface. For example, pillared superhydrophobic surfaces have been typically modeled using a simulation box containing a single unit cell consisting of the repeated motif composing the actual system: one pillar and the associated part of the bottom wall (see Fig. 1). hemeda2014; zhang2014; prakash2016; tretyakov2016; panter2016 This simulation setup prevents the breaking of translational symmetry of the liquid meniscus in the directions of the surface during wetting, which have been observed in a number of recent experiments. butt2013; lv2014b This limitation is especially important because symmetry breaking, and its origin, whether intrinsic of Cassie-Wenzel transitions giacomello2012b or due to the presence of dirt on the surface, DuanReview2016 is debated in the literature.
The aim of this work is to analyze in detail the effect of the choice of CVs and the simulation box size on the mechanism and energetics of collapse of superhydrophobicity in 3D textured surfaces with interconnected cavities, namely on pillared surfaces. We focus on this type of surfaces because they are largely used in technological applications verho2012 and in experiments, checco2014; butt2013 are simpler than other 3D surfaces to investigate, and yet preserve the main aspects of a large class of surface textures, namely that the surface cavities are interconnected, which makes this research of general interest. Moreover, the present results could be of interest also for other fields in which the density field is relevant to describe thermally activated events: nucleation of crystals, Kelton:2010fs; page2006; page2009 protein folding, miller2007 membrane poration, fuhrmans2015; smirnova2015 etc. In these fields there is an interest in developing coarse-grained models, which reduce the computational cost while capturing the relevant physics. A more detailed analysis of the physics of wetting and dewetting on nanopillared surfaces is contained in a related article. amabiliPRF
Concerning the CVs, we consider the (coarse-grained) density field at various levels of coarsening. The reason to focus on this CVs is manifold; first of all it has been shown to be adequate to properly describe the wetting path in atomistic simulations of simpler systems. giacomello2015 Secondly, the density field is the main variable of the sharp interface model and of the (classical) density functional theory, two continuum theories which have been successful at describing wetting/dewetting in confined systems talanquer1996; Talanquer:2001; Lutsko:2008 (see also Refs Kelton:2010fs; Meloni:2016it). Finally, by changing the level of coarse graining of the density field one can explore the two opposite limits, the highly coarse-grained case of the average fluid density in the surface corrugations, equivalent to the simple CV used in Refs. giacomello2012; savoy2012a; amabili2016; prakash2016, and the case of the density field computed on a grid with spacing of few atomistic radii or for continuum liquids. ren2014; panter2016; pashos2016
To compute the wetting path, i.e., the sequence of coarse-grained density fields (the CVs set) that the system takes on along the Cassie-Wenzel transition, we use the zero temperature string method. 111In Ref. giacomello2015 we have discussed the importance of finite temperature effects in the wetting mechanism. In particular, we focused on different wetting paths and on the overlap of the related reaction tubes, i.e., the region of CV space where most (a prescribed fraction) of the reactive trajectories pass through. While this effect could be important also for pillared systems, the methodological aspects investigated in this work do not depend on this type of effects. In fact, their discussion would increase the complexity of the work without shedding new light on the choice of CVs discussed here. maragliano2006; weinan2002 The string method allows one also to determine the energetics along the transition path. In the case of a single CV, the average density of fluid particles in the corrugations, the string method reduces to a simpler techniques, the restrained MD method (RMD TAMD). As we will discuss more in detail below, results obtained by RMD (or, equivalently, by temperature accelerated MD – TAMD TAMD) are analogous to those that can be obtained with umbrella sampling (US TorrieValleau), including its indirect version (INDUS Patel:2010cw), and boxed molecular dynamics (BXD Glowacki:2009hqa): simulation artifacts are not due to the specific simulation method employed but to the choice of CVs. Thus, the conclusions we will draw on RMD have a broader scope: they do not refer to a specific simulation technique but to the use of suitable or inadequate CVs.
Concerning the system, we considered two cases: one comprising a surface with square pillars and one with a single pillar, i.e. only one unit cell (Fig. 1). In both cases periodic boundary conditions are applied in the direction of the surface, and . The second case is analogous to the one adopted in Refs. prakash2016; zhang2014; panter2016.
Summarizing, in this work we investigate possible artifacts arising from the choice of the CVs and from the size of the computational sample, and the correlations between these two simulation parameters. The systems investigated are shown in Fig. 1.
Understanding the effect of the coarse graining and system size on simulation results, whether or not they introduce any sizable artifact, is crucial to assess theories on the wetting and recovery mechanism and energetics and to provide a practical guideline for future simulations of this process on even more complex surface topographies. This will help moving forward toward the in silico design of superhydrophobic surfaces with tailored properties. Such properties indeed depend on the surface characteristics (topography and chemistry) via the wetting/dewetting path, which should therefore be captured accurately in simulations.
ii.1 Atomistic model of the composite system
The atomistic model used to investigate the Cassie-Wenzel transition consists of a Lennard-Jones (LJ) fluid in contact with a nanopillared LJ solid surface (see Fig. 1, left panel). The fluid 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 (see below for details). Fluid and solid particles interact via a modified LJ potential
where and are the characteristic length and energy of the potential, respectively. is the interatomic distance between particles and . is the parameter controlling the hydrophobicity cottin2003 of the solid material; here we use , to which corresponds a Young contact angle 222The Young contact angle is the angle formed between a flat surface and the tangent to a sessile liquid droplet deposited on it at the liquid/solid contact point. . giacomello2012; amabili2015
The fluid particles are kept at constant temperature via a Nosé-Hoover chain thermostat. martyna1992 The pressure of the liquid is controlled through the application of a constant force on the particles of the upper wall, which acts as a piston, while the bottom wall is kept fixed during the simulation. amabili2016 Within the quasi-static simulation approach used in the present work, the vapor pressure depends, essentially, only on . amabili2016Yeomans Thus, the choice of and determine the driving force of the wetting process, , which is in the present case.
The two investigated systems are shown in Fig. 1. One consists of a grid of pillars, and the other of a single pillar, denoted as in the following. Periodic boundary conditions are applied in the and directions. The pillars are tall and thick, and are at a distance of from each other.
ii.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 observables that are function of the positions of the fluid particles. From a macroscopic perspective, the Cassie-Wenzel transition can be described by the advancement of the liquid/vapor meniscus (sharp interface models) and its twisting and bending. lv2014b; butt2013 An alternative macroscopic approach, the (classical) density functional theory (DFT), evans1992; lowen2002 considers the fluid density field as the fundamental quantity to describe wetting. ren2014; pashos2015; panter2016 It is therefore rather natural to use the coarse-grained density field as CV in atomistic rare event calculations. In practice, the region of the simulation box comprising the pillars is divided into parallelepipedic cells, and for each cell one counts the number of fluid particles inside it:
where runs over the fluid particles, 333, the characteristic function, is equal to if is in -th cell and otherwise. can be expressed as the product of two Heaviside step functions per Cartesian direction: , with and the begin and end of the k-th cell in the direction, and 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, miller2007 convenient, but others are possible. is the (smoothed) characteristic function of the -th cell, and is the volume of the cell.
In a quasi-static representation of the wetting process a key quantity is the probability density that the collective variables take on the values . corresponds to a macroscopic state of the system, i.e., the Cassie, Wenzel or any intermediate state visited during the transition. is frequently expressed in the logarithmic form and in thermal energy units, the Landau free-energy :
where is the inverse thermal energy and is the Dirac delta function. The probability density function is expressed in terms of the microscopic measure . High-probability states, i.e. Cassie and Wenzel, are global or local minima of .
allows to gather information on the mechanism and kinetics of the process. 444Indeed, to obtain a complete understanding of the transition mechanism and kinetics one needs also the committor function , i.e. the probability that a trajectory in reaches the final state first. In the string method discussed in this section, which is based on the backward Kolmogorov equation, the committor function is taken into account implicitly.vanden2006transition In principle, one could run a long MD and compute 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. maragliano2006 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 fields of the fluid along the Cassie-Wenzel transition. The fluid can follow many wetting paths among which the one identified by the string method is the most probable one. In the presence of large barriers other paths significantly different from the one identified 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. laio2005 The string method is, therefore, well suited to deal with the large set of CVs employed in this work.
The most probable path , with the path parametrization , satisfies the condition
where indicates that Eq. 4 refers to the component of orthogonal to the path. is the gradient of the free energy and is the metric matrix of elements
is associated to the change of variables from to . maragliano2006 is the independent variable of the parametric representation of the wetting path . 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 , and 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 the effective force has zero components orthogonal to the path.
In practice, the path is discretized in images, , with a constant distance between neighboring images,
The most probable path is obtained using an iterative algorithm, the string method weinan2002 (see App. A), which, starting from a first guess, produces a path satisfying the discretized version of Eq. 4. The main ingredients of Eq. 4, and , can be computed using, for example, RMD (Sec II.3).
The first guess path is obtained from a high pressure wetting trajectory. The wetting barrier decreases with the pressure and, at this pressure, a wetting event has been observed over a timescale of timesteps. Along this spontaneous wetting trajectory one extracts a number of configurations ( in the present case) and computes the associated coarse-grained density field, which represent the first-guess path. This approach, in which the string method is initialized from a continuous atomistic trajectory, prevents potential problems in the calculation of the distance between neighboring images due to the translational symmetry of the system: for example, in the simple definition of the distance used here, two configurations corresponding to a translation of the density field by one or more period in the and/or direction would result into two different values of , even though they represent the same physical system. The correct choice between these two configurations is the one which is the push-forward in time of the previous one. Since, however, the string is a local optimization algorithm (App. A), the initialization procedure via a continuous dynamics described above leads automatically to the selection of the proper density configuration 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 along the path, thus obtaining . Usually the free energy along the wetting is reported as a function of the liquid fraction , with , , and the total number of particles in the textures (the region used to define 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. S4 in the Supporting Information), it is possible to report the free energy in the more convenient and illustrative parametrization .555We remark that the passage from the coarse-grained density field 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 field but the opposite –that there is only one value of the density field 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 field , and does not mean that and , i.e., CVs at different 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.
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 field in the present case) can affect 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. bolhuis2002; maragliano2006 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 . The transition interface is locally defined 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 fluid toward the Wenzel or Cassie state from the putative transition state, the committor test has been rarely performed. giacomello2015 Here we introduce a simpler test allowing to check whether the CV set satisfies minimal necessary conditions, namely that 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.
ii.3 Restrained MD and related simulation techniques
In order to perform one iteration of the string algorithm one needs to compute and at the set of points . In this section we show how these quantities can be computed by RMD.
Consider Eq. II.2 defining the Landau free energy and replace the Dirac delta functions on the r.h.s. by smooth approximations of Gaussian form . Within this approximation each component of the gradient of the free energy reads
Thus, is approximated by the expectation value of the observable over the conditional probability density function . Here, we have set , which is a trade-off between the convergence of 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 is the measure of the canonical ensemble, can be cast in the form , which suggests that the conditional probability density can be sampled by a constant temperature MD driven by the augmented potential , the so-called Restrained MD. Indeed, this argument can be extended to the isothermal-isobaric ensemble, provided that one uses a molecular dynamics suitable to sample this ensemble.
In practice, is computed as the time average of along the RMD. can also be computed as a time average along the same RMD, this time considering the observable .
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 final state. Thus, in this case RMD allows to compute directly the free energy via the numerical integration of . In this sense, sometimes we refer to RMD as a proxy for a single-CV string.
Iii Results and discussion
This section is divided in three subsections, i) pillars system (Sec. III.1) and ii) pillar system (Sec. III.2), in which we report and discuss results of the effect of different levels of coarse graining of the density field on the wetting path and its energetics for each sample independently, and iii) Size effects: comparison between and pillars systems (Sec. III.3), in which we discuss the effect of the size of the surface sample on the results.
iii.1 pillars system
Figure 3 shows the evolution of the liquid meniscus along the collapse of the Cassie state for ( mesh with a cell of ), ( mesh with cells of ), and ( mesh with cells of - 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 identified using the Gibbs dividing surface, i.e., the locus of points where the fluid density coincides with the mean value between the liquid and vapor bulk values, . Thus, identifying the meniscus from an ensemble of atomistic configurations requires computing a coarse-grained density field, . For this computationally inexpensive post-processing operation we use a finer level of coarse graining than that used for the CVs, namely a points grid (cell ).
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 fill the inter-pillar space, with the meniscus assuming an essentially flat 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 in Fig. 5(a).
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 flat meniscus the liquid-solid surface increases resulting in an increase of the energetically unfavorable solid-liquid contribution; ii) when the liquid starts spreading over the bottom wall the two liquid-vapor and vapor-solid interfaces are replaced by a single liquid-solid one, which results in a reduction of the free energy. Thus, the free energy switches from an increasing to a decreasing trend when the meniscus comes in contact with the bottom wall.
Step B of the wetting path is significantly different in the three cases (Fig. 3). With CV, the liquid touches the bottom wall by completely filling the space between two rows of pillars. Due to the 2D periodicity, this amounts to having an infinite strip of liquid touching the bottom of the textured surface parallelled by an infinite strip of vapor, with the corresponding meniscus facing the bottom wall at a distance of ca. . The process proceeds with the liquid progressively invading neighboring rows or squares of pillars.
For both and 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. However, wetting evolves in a rather different way. In the case of CVs the liquid percolates in neighboring inter-pillar regions; with CVs it wets disjoint regions.
A remarkable aspect of step B with and CVs is that both present discontinuities in the path corresponding to abrupt changes of the meniscus configuration (Fig. 3). For example, in the case of CV the filling of the corrugations proceeds with the liquid suddenly wetting entire squares or rows among pillars. Concerning the 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 and CVs to distinguish among different meniscus configurations: more than one configuration of the meniscus is possible for a given value of the CVs.
Step C of the path is similar for all the cases: a rarefied 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.
Concerning step B of the wetting path, which includes the kinetic bottleneck of the process – the transition state, an important question is whether the paths identified with the three sets of CVs are all physically sound. Before addressing this question let us 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. maragliano2006 Thus, for systems evolving according to continuous dynamics, one would expect that observables change smoothly between successive images of the discretized path. Of course, in all collective-variable-based rare event techniques, the CVs used to follow or accelerate the process change by a prescribed value between images; for example, in the present string calculations the distance in the value of CVs between successive images is constant by construction, Eq. 6. However, one expects that this property holds for all the fundamental variables of the system. In particular, in the case of wetting of textured surfaces, which is well described by the (classical) density functional theory, ren2014; tretyakov2016 a necessary condition for the CVs to be suitable to describe the process is that the density field (computed, virtually, on an infinitely resolved grid) changes in a continuous way along the path regardless of the degree of coarse graining considered. To verify whether this condition is met, we computed the nondimensional Euclidean distance between (ensemble averaged) density fields computed on a very fine mesh (cell ) at successive images along the paths identified via the three sets of CVs:
where the normalization is performed with respect to the bulk liquid density , and the sum runs over the cells indices of the very fine mesh, which is the same for all the CVs. We remark that, apart for the normalization and the mesh used in Eq. 8, is equivalent to the normalized arc length (Eq. 6). The value of far from discontinuities depends on the absolute arc length of the path in the CV space and on the number of images in which it has been discretized. Thus, there is no trivial relation between the number of CVs and the value of . In practice, for the case of the wetting of pillared surfaces, far from discontinuities takes similar values for all the considered sets.
It is observed (Fig. 4) that in steps A and C, independently of the number of CVs, indicating that there are no major differences in the continuity of the paths identified by the three CVs in these regimes. On the contrary, in step B the paths obtained with and CVs show discontinuities with a sizeable increase of , which reaches and in the first and second case, respectively. These results suggest that and CVs are insufficient 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.
The difference in the wetting paths is reflected in the difference among the free energy profiles (Fig. 5(a)). In the case of CV we observe an extended, relatively flat domain of the free energy profile in step B. Indeed, careful analysis shows that this region is composed of two linear segments with a different slope (Fig. 5(a)). The discontinuity of the first derivative is observed in correspondence of the transitions from one morphology to another (see inset of Fig. 5(b)), namely i) from the flat meniscus to the liquid partly wetting the bottom wall (transition from step A to B); ii) from the liquid completely wetting one row of pillars to two rows; 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 profile is more regular in the case of and CVs, with a well defined maximum of the free energy between the Cassie and Wenzel states.
The quantitative comparison of the three free-energy profiles reveals two important aspects. First, the difference of free energy between the Cassie and Wenzel states, , which determines their relative stability, depends on the CVs. and CVs yield consistently , while the CV case is ca. higher. This discrepancy, in turn, affects the barrier of the Wenzel-Cassie transition, reducing its estimate by a corresponding amount in the case of CV. Second, also the Cassie-Wenzel free-energy barrier changes from one CVs set to another. In particular, is approximatively and lower for the and CVs, respectively, as compared to the CVs case.
We remark that the differences in the free energy profiles with the three CV sets is not due to an insufficient discretization of the RMD and string paths. To show this, we repeated the RMD calculation in B and C with three times the number of points (see Supporting Information Fig. S5(b)). The original and more accurate profiles show minor quantitative differences. To further confirm that the difference in the free energy profiles does not depend on an insufficient discretization of the path but on the difference of the CV sets employed, in the region across one of the changes of slope of for 1 CV we increased the number of point by one order of magnitude (Fig. S5(a)). Also in this case we observed minor quantitative differences in both the free energy profile and mean force.
iii.1.1 Dependence of on the CVs
The theoretical arguments developed in App. B indicate that, the values of should be close for the three CVs; the present results conflict with this theoretical prediction. The origin of this mismatch can be illustrated considering a model 2D polynomial potential which has an analytical expression and on which different CVs can be analyzed in detail (Fig. 6(a)). In panel b) the profile of 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 profile obtained with higher coarse graining; this 1D potential is obtained by projecting along the axis according to Eq. II.2: , where is the marginal probability density function.666The marginal is the probability density function that takes a given value while takes any value. Of course the concept of marginal is not limited to the case of two random variables but can be extended to the case of any number of random variables, and refers to the event that a subset takes a given value while the rest of the variables take any value. Figure 6(b) shows that , the energy difference between the two minima, is indistinguishable if computed along the string path or using the 1D potential ; the differences along the path will be analyzed later on when the barriers are considered.
We further compute as one would obtain from an RMD simulation (thought RMD - tRMD) performed using as CV; this is expected to provide a numerical approximation to . In tRMD one imagines to start from the minimum increasing step by step . In finite time, for each , the system explores configurations within few from the local minimum at the constrained value of belonging to the basin of (green curve in Fig. 6(a)). When the barrier along separating the and basins is sufficiently small the system jumps into the other basin and successively oscillates around the other branch of the green curve which brings to the basin . In Fig. 6(a) the transition from the initial to the final branch of the path occurs when the barrier along the direction is . From the ensemble of configurations available at each one estimates the mean gradient
from which the free-energy profile is reconstructed by integrating with respect to .
As noted, the values of obtained from and the projected coincide, while is not captured correctly by tRMD. We remark that the relative error on estimated via tRMD is , of the same order of magnitude of the error on with 1 CV. The reason of the error on is that the mean gradient used to compute the free-energy profile in tRMD, Eq. 9, differs from the actual one near the jump between the two valleys of the potential (see the inset of Fig. 6(b); additional details are discussed in App. B). In particular, 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 difference in the between and or CVs is due to the simulation protocol. The transition from one basin to another in the 2D potential is equivalent to the meniscus jumps observed in the wetting process with CV. Thus, we believe that the difference in between and or CVs can be ascribed to the severe discontinuities in the gradient of the free energy in the former case, which are induced by the extreme coarse graining of the associated density field in the case of one single CV. These discontinuities are shown for the case of a single collective variable in the inset of Fig. 5(b).
Indeed, this problem, which is also at the origin of the hysteresis observed with other rare event techniques when one uses insufficient or inadequate CVs, could be solved combining RMD with parallel tempering, 1986PhRvL..57.2607S which allows to correctly sample the conditional probability in the space complementary to the collective variables (see, e.g., Refs. Orlandini:JournalOfStatisticalPhysics:2011 and PhysRevB.83.235303 for the application of this combined approach to phase transitions). This approach would, thus, bring of CV closer to the value computed with or CVs; however, to the best of our knowledge, due to the high computational cost this approach has been rarely adopted in the context of wetting transition. giacomello2012
iii.1.2 Dependence of the wetting barrier on the CVs
In the following, when comparing the free energy barrier of different CVs, we will first focus on the zero temperature limit and we will then consider finite temperature effects. It must be stressed that zero and finite temperature refers to the CVs degrees of freedom, while the atoms, and thus the free-energy landscape, are always at the physical temperature . CV is a subset of CVs which, in turn, is subset the CVs. Thus, in the present context, the zero-temperature limit of the CV set means that all the remaining uncontrolled degrees of freedom to complete the -dimensional CV space take the value corresponding to the constrained local or global minimum at the present value of the restrained CV. An analogous argument holds for the comparison of the barriers with and CVs, and vs CVs. The reason to split the comparison in the zero temperature limit and finite temperature effects is that this allows to carry on an accurate theoretical analysis, bringing to some interesting conclusion on the possibility or impossibility to establish an ordering of the barriers with the size of the CV set.
Let us consider first 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 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 (blue line). The string in variables departs from the zero-temperature tRMD path in the intermediate region, undergoing the transition from the initial to the final basin at a different point. Using a different language, the trajectories with one and two variables cross the separatrix, the surface separating the two attractive basins, at different points. In Ref. maragliano2006 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)). We remark that this ordering of the barriers is valid only in the zero-temperature limit, and is due to the fact that the system keeps moving along the line of constrained minima of the basin beyond the point at which the potential of the minimum in the basin would be lower.
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 App. B). One notices that the barriers with and CVs are lower than that with , which conflicts with the zero-temperature analysis above and the intuitive expectation that the barriers should decrease as the number of CVs increases. This suggests that there must be other effects which are not present in the zero-temperature limit.
At finite temperature there are two additional effects to take into account. The first concerns the point where the transition occurs with a reduced number of CVs. At finite temperature this is attained where the orthogonal barrier is of the order of several , 777The transition from the initial to the final branch of the path is performed when the barrier along the direction is zero and for tRMD and tRMD, respectively. The first condition is met at the green arrow of Fig. 6(a). The second condition occurs very close to this point. i.e., when the transition time is of the same order of magnitude of the string/RMD simulation time. This implies that the transition occurs before than in the zero temperature case, i.e., closer to the initial state, which results in a reduction of the barrier.
The second effect 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 configuration space the system can explore and, thus, the lower should be the configurational entropy at a given point along the path. For example, in the case of 864 CVs at each point along the string the fluctuations allowed the density field to vary over a spatial scale smaller than the cell in which the system is partitioned, i.e. on a spatial scale of . In the case of CVs the density can fluctuate on a much larger spatial scale of . This effect can be better illustrated in the case of the 2D polynomial potential. At each point along the string in the system can take only one configuration, and the entropy is zero. On the contrary, in the case of finite temperature tRMD the system oscillates in the direction around the minimum at the current value of ; the entropy at each point along the RMD path is different from zero.
This effect is qualitatively valid in general, at all points along the path. However, the value of the entropy depends on the profile of the potential in the space orthogonal to the CVs at the current point along the string/RMD. Thus, in the 2D potential the entropy of the system at a given depends on the shape of in the direction: if the potential is stiff the entropy is low, if it is shallow the entropy is high. For example, the potential in the direction is stiff at the initial and final states and shallow at the transition state: thus the entropic contribution, , tends to reduce the barrier as compared to tRMD at (see inset of Fig. 6(b)).
The combined effect of the temperature in reducing the barrier, anticipating the transition and by entropic effect, is illustrated by comparing tRMD at zero and finite temperature. 888It is worth remarking that there is no obvious ordering among the free energy profiles of the single variable because the value of the free energy depends on several competing factors: the value of the free energy at the constrained minima in the two basins, the associated curvature and the temperature (see also App. B). Indeed, Fig. 6 shows that at finite temperature the tRMD barrier can be lower than the string one.
A similar scenario is likely to occur for the wetting process, in which and CVs might have a higher entropy than at the transition state. For the case of CVs this hypothesis is supported by the visual inspection of the trajectories. The fluctuations of the meniscus close to the Cassie state, in which the liquid/vapor interface is pinned at the pillars top, is small in all cases. On the contrary, at the transition state the meniscus shows larger fluctuations for CVs than for .
To validate this hypothesis we have computed the wetting entropy by subtracting from the free energy the enthalpic term, given by the average Hamiltonian at each image of the path plus the pressure term (, with denoting the ensemble average of the Hamiltonian at the present point of string/RMD and and the pressure and volume, respectively, of the phase , liquid or vapor). 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 figure). In correspondence of the transition states both the and CV cases present an increase of the entropy, which is, however, more pronounced for the CVs. This, indeed, explains why the CVs set has a barrier lower than the one.
The scenario is more complex for the 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 fluctuations, contributes to the entropy of the system. This is different from the other two cases, in which the liquid-vapor interface is eliminated gradually, and the curved liquid-vapor interface joining the triple line at the bottom wall with the flat part of the meniscus increases the entropy. At the boundary between B and C, when teh CV system changes from a well defined liquid-vapor biphasic system to one made of a liquid plus an extended and diffused liquid-solid interface, the entropy suddenly increases. The final part of the entropy profile is similar for all cases: it decreases following the absorption of the regions of rarefied liquid (Fig. 3).
Summarizing, the level of coarse graining affects the wetting barrier in three different ways: i) at zero temperature a lower degree of coarse graining increases the barrier because the separatrix is crossed outside of its minimum; at finite temperature ii) the crossing of the separatrix is anticipated and the barrier reduced in all cases, but the effect is higher for more coarse grained density fields and iii) entropic effects may further reduce the barrier for more coarse grained density fields. However, this last argument has exceptions for too coarse grained density fields, e.g., a single CV, which changes the wetting mechanism toward less entropic paths. The balance between these three effects cannot be easily predicted. Rather, they simply explain the ordering of barriers observed with , and CVs.
Let us close this section drawing some conclusions on the effect of the degree of coarse graining on the study of the wetting/dewetting process. Extreme coarse-graining also leads to severe errors in the estimation of , up to ; in addition, it also causes discontinuities in the wetting path. Extreme coarse-graining affects the wetting barrier, which is underestimated by as compared to the finer density field case. Overall, this suggests that quantitative predictions require a careful choice of the CVs. We remark that these effects 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 App. C.
iii.2 pillar system
For the pillar system we perform an analysis similar to the one of Sec. III.1. In this case we considered and CVs, having the same grid spacing of the and CVs cases of the pillars, respectively. Fig. 8 reports , , and . The general characteristics of the wetting mechanism observed for the pillars system are preserved also in the smaller sample: the liquid initially depins from the pillars top, slides along their side with a flat meniscus (A), then touches the bottom wall in a point between a pair of pillars (B), and finally the vapor bubble is absorbed (C). It is worth remarking that analogous mechanisms have been reported in the literature for the case of water wetting a pillared surface. In this case, simulations were performed using INDUS using 1 CV analogous to the one used in the present work. prakash2016 The CVs case presents a smooth transition from one regime to another, with remaining almost constant at . On the contrary, with 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 and CVs free energy profiles present no qualitative differences: both curves are characterized by a single, well defined maximum in correspondence of the configuration in which the meniscus touches the bottom wall. Concerning quantitative aspects, while is the same in both cases, the barrier of CV is lower than for CVs. This is due to the entropic effect explained in the previous section.
A more careful analysis of the free energy profile close to the Wenzel state reveals differences in the curvature of the free energy as a function of the filling fraction. In order to study such differences and to obtain a more accurate free energy profile in region C, the path with CVs in the range has been discretized in images (Fig. 8(c)). With CV a parabolic profile extending over a broad range of filling fraction is observed, which is generally attributed to Gaussian fluctuations of the liquid density as discussed in the Lum-Chandler-Weeks theory of hydrophobicity. lum1999 A similar trend has been observed by other authors using similar CVs . remsing2015; amabili2016; prakash2016; giacomello2012 On the contrary, with CVs the parabolic trend is observed over a much narrower domain, after which the free energy profile shows a dependence of the type typical of macroscopic theories, which accounts for the liquid-vapor and solid-vapor interface energies giacomello2012a (see Supporting Information for further details). This is reflected on an earlier formation of the liquid-vapor meniscus, as shown in the snapshots of Fig. 8(a). The wetting mechanism identified with CVs is characterized by a strong initial () depletion of the density of the liquid in the corrugation before the meniscus is formed. On the contrary, the path obtained with CVs indicates a more modest density depletion before the liquid-vapor interface forms.999We remark that similar results are present also in the pillar case, in which a different slope of in the domain C is apparent (Fig. 5). The detailed analysis of the early stage of dewetting, which requires a fine discretization of the path with a large number of string images, was performed for the case which is computationally more convenient. Considering that the CVs is a superset of CV (overall density), i.e., the CV-space includes and extends the CV one, the different behavior of the free energy indicates that the overall density is inadequate for describing the wetting and dewetting path. Since differences arise already at the beginning of dewetting, relatively close to the Wenzel state, one must conclude that also in the presence of moderate dewetting barriers, which occur closer to the fully wet state, the prediction of paths and rates obtained from simulations based on a single CV might be inaccurate.
In several works remsing2015; prakash2016; amabili2016 the use of the overall density as the only descriptor of the wetting/dewetting process has been justified on the basis of the Lum-Chandler-Weeks theory lum1999 of hydrophobicity. However, the above observation that CV is insufficient to characterize thermally activated wetting/dewetting processes suggests that the extension of the theory of the hydrophobic effect out of its original equilibrium scope may bring to overlooking important aspects of the transition. Indeed, fluctuations 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.
iii.3 Size effect: comparison between and pillars systems
In the previous sections we focused on the effect of the choice of CVs on the wetting mechanism, energetics and kinetics; here we concentrate on the effect of the size of the surface sample, the number of pillars in the and directions. To this end we compare results obtained with CVs for the surface and CVs for the one. At the transition state the two systems present qualitatively similar configurations. In both cases the liquid touches the bottom wall at a point between two pillars (Fig. 9). However, the surface cannot reveal the complex nature of the final part of the wetting. Indeed, for the case we observe the formation of a percolating (random) network of vapor bubbles, which is absorbed when the system approaches the Wenzel state. amabiliPRF At variance with this scenario, for the 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, and (compare top and bottom panels of Fig. 9).
Concerning the quantitative comparison of results, Fig. 10 reports the free-energy profiles for the and pillars systems. Given the different size of the two systems, such an analysis is perfomed multiplying by a factor the free-energy profile of Fig.8b. While, as expected, the free-energy difference between Cassie and Wenzel states, , are in agreement, the barrier of the pillar system is higher than the one. Moreover, also the position of the transition states is different in the two cases: for the pillars system as compared to for the one. This is the result of small size effect, more specificallt, of the (periodic) repetition of the bent meniscus in the replicated system, see Fig. 9. pillars system. A fairer comparison can be obtained considering the free energy of the transition state plus eight times that of the last flat meniscus configuration of the system. However, also in this case the barrier is 180 higher than the . This brings us to the conclusion that the system is affected by important finite-size errors, due to the enforcement of an artificial symmetry.
In this work the Cassie-Wenzel transition has been studied using rare event methods, namely the string and RMD techniques. The collective variable used is the density field 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 field, often used in the recent literature, introduces qualitative and quantitative artifacts. These artifacts include the identification of a sequence of physically disconnected density fields for the Cassie-Wenzel transition path as well as wrong estimates of the relative energy between the Cassie and Wenzel states and the barrier separating them.
We have also investigated the effect 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 ( pillar system) is insufficient when dealing with interconnected textures. Indeed, the periodic boundary conditions impose a non-physical translational symmetry of the meniscus leading to a 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 finger 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 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.prakash2016 However, the detailed physics of the wetting/dewetting processes could be qualitatively different from the simulated one and the value of the free-energy barriers could be off by tens of . 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 different cases in order to devise design strategies, e.g., for preventing the complete wetting or helping the recovery of the Cassie state. A posteriori, one could evaluate the choice of CVs: typical symptoms of the deficiency of the adopted description is the presence of non-differentiable points in the free-energy profiles or the appearance of jumps in the density along the wetting/dewetting processes.
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. . We acknowledge PRACE for awarding us access to resource FERMI based in Italy at Casalecchio di Reno.
Appendix A The string algorithm
The string algorithm aims at finding a path in the CV space satisfying the condition that the force orthogonal to the path, via a metric matrix , is zero (Eq. 4):
In practice, the path is discretized in images, and takes values with . As already mentioned in the main text, the images are enforced to be at a constant distance from each other (Eq. 6).
One can introduce the following artificial dynamics
where is an artificial time bearing no direct physical meaning. For , converges to the solution of Eq. 10. The r.h.s. of Eq. 11 can be computed projecting out the component of the force parallel to the path at the present :
with tangent to the path and its transpose. The numerical determination of the tangent to a discrete path can be inaccurate, which brings to numerical instabilities in the iterative procedure of Eq. 11. Instead of removing the tangential component of the force, is first evolved according to the entire force, then the equidistance among images is restored. e2007 This second step is performed by moving the images along the polygonal curve interpolating the path. This renormalization step, by removing the effect of the tangential component of the force, is formally equivalent to applying the projection.
Appendix B Relation between the probability density function of different collective variables
Direct comparison of the free energies of different sets of collective variables is, in general, not possible. However, in the case of the sets of CVs used in this work, and because of the steep profile of the free energy, this comparison is, to a large extent, meaningful. The reason of this is explained in the following.
Consider first the cases of and variables, the associated free energy and and probability density function and . To make the argument more straightforward, assume that the relevant ensemble is the canonical one; in this case
At low temperature, if the free energy has a minimum at each value of , , the integrals in Eq. 13 can be solved to a good approximation via the Laplace’s method. This amounts to replacing the integrand with the exponential of the second order expansion of around . Thus, , with the second derivative of the free energy with respect to at the minimum. The corresponding formula for the free energy reads .
The free energy is typically written with respect to a reference state, say the state associated to the global minimum of the free energy , and the relative free energy is . Henceforth the relative free energy of the state is the same as the one of the two-dimensional free energy at plus a term which is proportional to the logarithm of the second derivative of the free energy. This second term is expected to be exponentially negligible as compared to the first and goes to zero for .
It is worth remarking that the above analysis holds true also in the case in which presents more than one minimum in at a given value of . In this case, one could divide the integral at the numerator of Eq. 13 in the sum of several integrals centered around each minimum, and each integral can be solved with the Laplace’s method. The results is . Given the exponential dependence of the probability density function on the free energy, , where is the absolute minimum at . When two or more minima have the same probability,
In order for the above analysis to be useful in the present work, it must be extended to the case in which the variable of which one wants to compute the marginal probability is a function of the variables of the original set, not just one element of the set. For example, the CV variable set discussed in the main text, consisting of the number of particles in the control volume, can be written as the sum of the collective variables of the largest CVs set. Consider first the transformation of variables , with . The Jacobian of this transformation is , which implies that . One can then use the argument of the previous paragraphs to show that , where is the configuration identified by the string/RMD calculation at .
The above arguments prove that at low the free energy profiles with all the CVs sets should be similar (provided that they are computed along the same path).
Appendix C Comparison among string/restrained MD, umbrella sampling and boxed molecular dynamics
In this work we have analyzed the effect of the choice of the level of coarse graining of the discretized density field on the mechanism and energetics of the wetting process performing string and RMD calculations. However, present results are not limited to the special rare-event simulation technique adopted: analogous artifacts and drawbacks affect simulations performed using other techniques with the same CVs. The most used techniques for the study of wetting have been umbrella sampling (US)TorrieValleau and boxed MD (BXD);Glowacki:2009hqa in this section we perform a comparative analysis of these methods and explain why one should expect analogous problems with these techniques when using the average density in a given region as CV.
US relie on an algorithm which is similar to RMD: they both consist in performing MD with an augmented potential of the form of . However, in US the parameter is, typically, set to much smaller values than in RMD. Thus, the system is allowed to explore a larger space as compared to RMD, in which the system visits only configurations consistent with the condition . In US one first computes the free energy within the range of collective variable values explored at the present value of . Then, the branches of the free energy profiles of different umbrella potentials are joined by imposing that they take the same value over the overlapping regions of CVs. More elaborate and efficient techniques exist to compute the free energy from US simulations WHAM but their discussion is beyond of the scope of the present article. In the presence of morphological transitions, when a value of the CV in one or the other basin corresponds to different configurations, the hypothesis of same probability density for same value of the CV is not satisfied, which might bring to large errors in the reconstructed free-energy profile (Fig. 11(a)).
In BXD there is no augmented potential, and the atoms evolve according to the physical potential . However, the system is confined within a box in the (single) CV space, defined by a lower and upper value of the collective variable. This box is equipped with reflecting boundary conditions, which allows to extensively sample regions of overall low probability. The calculation of the free energy profiles in BXD is a three-steps procedure; first, one determines the free energy difference between neighboring boxes from the rate of transition of forward and backward trajectories between them: , with and rates of transition from box to box and vice versa, respectively. From the complete set of for all pairs of consecutive boxes one can compute the probability that the system is within a given box. Secondly, the probability density of the CV within each box, , is computed from the histogram of the corresponding trajectory. Thirdly, the Landau free energy is obtained from the logarithm of the total probability , i.e., via Eq. II.2. When the system is characterized by two disjoint initial and final basins as in Fig. 11(b), the flux of trajectories across a given interface is different if the system is in one or in the other. Indeed, also depends on which attractive basin the system is in. Thus, also in this case, inadequate CVs may yield wrong free energy profiles.
Summarizing, the error in the estimation of the free energy is not associated to a specific rare-event technique (RMD, US, BXD, etc.), rather, it depends on the choice of CVs.