{'fraction_non_alphanumeric': 0.03781085737002518, 'fraction_numerical': 0.014201157294933521, 'mean_word_length': 4.446114986769305, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Characterizing Surface Wetting and Interfacial Properties using Enhanced Sampling (SWIPES)

We introduce an accurate and efficient method for characterizing surface wetting and interfacial properties, such as the contact angle made by a liquid droplet on a solid surface, and the vapor-liquid surface tension of a fluid. The method makes use of molecular simulations in conjunction with the indirect umbrella sampling technique to systematically wet the surface and estimate the corresponding free energy. To illustrate the method, we study the wetting of a family of Lennard-Jones surfaces by water. We estimate contact angles for surfaces with a wide range of attractions for water by using our method and also by using droplet shapes. Notably, as surface-water attractions are increased, our method is able to capture the transition from partial to complete wetting. Finally, the method is straightforward to implement and computationally efficient, providing accurate contact angle estimates in roughly 5 nanoseconds of simulation time.

liquid over its vapor, can be characterized by estimating either k or θ .
The most common approach for characterizing wettability is the direct measurement of contact angle using the socalled sessile droplet method, wherein θ is determined from the geometry of a liquid droplet supported by a solid surface. Although the sessile droplet method is usually associated with the experimental determination of θ , molecular simulations have also made extensive use of this method. Indeed, both molecular dynamics (MD) and Monte Carlo (MC) simulations have been used to estimate contact angles for diverse interfacial systems, ranging from toy models [6][7][8][9][10] to realistic molecular models such as water/minerals [11][12][13][14][15][16][17] and textured surfaces [18][19][20] . However, in contrast with the millimeter-sized sessile droplets used in experiments, the droplets used in molecular simulations are nanoscopic. At the nanoscale, the geometry of spherical droplets can be influenced by the tension associated with the three-phase contact line. As a result, estimates of θ can depend strongly on the size of the simulated droplet 21,22 . Such line tension effects can be mitigated by using a cylindrical droplet that is infinitely long (due to periodic boundaries), because the three-phase contact line then has a fixed length and is independent of droplet size 23 . However, a number of studies have recently shown that θ extracted from the geometry of such cylindrical droplets nevertheless depends on the curvature of the droplet [24][25][26] . In addition to the challenges posed by line tension and finite size effects, the determination of θ from droplet geometry is also plagued by a certain degree of ambiguity regarding the exact location of the solid-fluid interface [27][28][29] .
Given the limitations of the sessile droplet method, a number of free energy methods have been proposed for estimating k. The interface potential approach involves estimating the free energy per unit area needed: (i) to wet a surface in contact with vapor, γ VL + γ SL − γ SV = γ VL (1 − k) (spreading potential), and (ii) to dewet a surface in contact with liquid, γ VL + γ SV − γ SL = γ VL (1 + k) (drying potential), thereby providing estimates of both γ VL and k 30 . Errington and co-workers have illustrated the utility of this approach, and have employed it extensively in conjunction with grand canonical Monte Carlo simulations to study the wetting properties of diverse surfaces by a number of fluids [30][31][32] . However, because using grand canonical Monte Carlo to simulate large or complex fluid molecules can be inefficient, and performing molecular dynamics simulations in the grand canonical ensemble is cumbersome, the interface potential method has not been applied to the study of such complex fluids. To this end, approaches that employ molecular dynamics simulations have been developed, e.g., the phantom-wall and dry-surface methods by Leroy and Muller-Plathe [33][34][35] . These approaches are similar in spirit to the drying potential described above, i.e., they employ processes that enable estimation of the work of adhesion. In the phantom-wall method, a structureless repulsive wall is employed to push the liquid away from the surface of interest 33,34 , whereas in the dry-surface method, the surfacewater attractions are reversibly turned off 35 . These methods have been widely applied to calculate wetting properties for a number of different systems 33,[36][37][38] . However, in each of the above free energy methods, the solid-liquid interface is replaced by both a solid-vapor and a vaporliquid interface 30,33,34,39 . Such formation of vapor-liquid interfaces in the presence of periodic boundary conditions involves transitions between dewetted morphologies [40][41][42] , which can give rise to hysteresis and complicate estimation of the associated free energies [43][44][45] .
Moreover, in the phantom-wall method, a solid-liquid interface is replaced not by a solid-vapor interface, but by a solid-vacuum interface. The corresponding free energy per unit area, γ VL + γ S,Vac − γ SL approximates γ VL + γ SV − γ SL = γ VL (1 + k) if the solid-vacuum surface tension, γ S,Vac ≈ γ SV . Such an assumption is reasonable when solid-fluid attractions are relatively weak, and there little to no adsorption from the vapor phase onto the solid surface. However, γ S,Vac is expected to differ substantially from γ SV when the solid-fluid attractions are significant. Indeed, Kanduc and Netz 46,47 have recently shown that for a sufficiently solvophilic surface, there is substantial adsorption from the vapor phase onto the solid surface, which lowers γ SV (relative to γ S,Vac ). To do so, the authors introduced a free energy method that enables estimation of the spreading potential, and is able to capture the transition from partial to complete wetting as solid-fluid attractions are increased. However, the method is computationally demanding and requires roughly 50 simulations for each estimate of k.
Here, we introduce a versatile and computationally efficient method, which we call "Surface Wetting and Interfacial Properties using Enhanced Sampling" or SWIPES. SWIPES makes use of the indirect umbrella sampling (IN-DUS) technique 48,49 to systematically wet the surface of interest, and estimate the free energy change associated with the process. It can be used in conjunction with molecular dynamics simulations to obtain estimates of not just wetting coefficients k, but also of vapor-liquid surface tensions, γ VL . We illustrate SWIPES by using it to study the wetting of Lennard-Jones (LJ) surfaces by SPC/E water 50 . We estimate wetting coefficients for surfaces with a wide range of solid-water attractions, and compare our results with those obtained from droplet geometries. We also show that our method is capable of capturing the transition from partial to complete wetting as the solid-water affinity is increased.
The rest of the paper is organized as follows: In section B, our approach for estimating wetting coefficients is described. Details pertaining to the molecular models used and the simulation methods employed are given in section C. In section D, contact angles calculated from free energies and from droplet geometries are compared, and the transition from partial to complete wetting as well as the computational efficiency of the method are discussed. Finally, we summarize our findings in the conclusions section.

B. Estimating Wetting Coefficients
In this section, we describe our approach for calculating the wetting coefficient, k, defined as: The wetting coefficient is related to the contact angle (θ ) through k = cos θ .

Outline of Our Approach
To estimate k for the wetting of a flat solid surface by a liquid in equilibrium with its vapor, we employ the schematic shown in Figure 1. As the number of waters, N v , in the observation volume, v, increases, and the system moves from state A to state B, the area of the solid wetted by the liquid increases by ∆A = 2L∆H, where L is the length of the simulation box in the direction perpendicular to the page, and ∆H is the distance that the vapor-liquid interface advances along the solid. The work, ∆F, that must be done to wet the solid as the system moves from state A to state B, is given by Although a (curved) vapor-liquid interface and three-phase contact lines are present in our setup, the corresponding interfacial area, its curvature, and the contact line length are all unaffected by the advance of the liquid slab along the solid surface, and thereby do not influence ∆F. The wetting coefficient (k) or the contact angle (θ ) can thus be estimated from the ratio of ∆F and ∆H as We note that the ratio ∆F/∆H ought to be the same (−2γ VL Lk), regardless the positions of vapor-liquid interface in state A and state B, as long as the solid surface under consideration is homogeneous, i.e., free of chemical heterogeneity or roughness.

Indirectly SamplingÑ v
To efficiently estimate ∆F/∆H, we employ enhanced sampling simulations that pull waters into v and wet the surface. In principle, this can be accomplished by biasing the number of waters, N v , in the rectangular observation volume, v, shown by the dashed lines in Figure 1. In practice, because N v is a discrete function of particle positions, R, biasing it using molecular dynamics simulations would result in impulsive forces. We therefore bias N v indirectly by instead biasing a coarse-grained number of fluid molecules, N v , following the indirect umbrella sampling (INDUS) prescription 48,49 ;Ñ v is closely related to N v , but is a continuous function of R. The precise definition ofÑ v and its dependence on R can be found in the ref. 49. Our choice of the location of v is discussed in further detail in the Supplementary Information.

Biased Simulations to Wet the Surface
We choose a biasing potential that is parabolic inÑ v (R), with κ and N * parametrizing the potential. We further choose state A and state B to correspond to biased ensembles with different values of N * , say N * A and N * B . Although this choice is convenient, it is certainly not unique; alternative definitions can be chosen for the two states, e.g., constant-Ñ v ensembles with differentÑ v -values. Nevertheless, this choice enables us to reversibly wet the surface, and to estimate the corresponding free energetics, which we characterize using the free energy difference between the biased and unbiased ensembles, F κ,N * . Similarly, to characterize the position of the vapor-liquid interface, H, we employ the biased ensemble average of the center of mass of water in the simulation box in the direction (x) perpendicular to the vapor-liquid interface, i.e., H ≡ x COM κ,N * , where O(R) κ,N * represents the average of O(R) in the biased en-semble. As shown in the Supplementary Information, both the position of the vapor-liquid interface, H, and the corresponding free energy, F κ,N * , vary linearly with N * ; i.e., both f ≡ dF κ,N * /dN * and h ≡ dH/dN * are independent of N * . Briefly, this stems from the fact that the balance of forces at the 3-phase contact line is governed entirely by the corresponding interfacial tensions, and the shape of the vapor-liquid interface is independent of the biasing potential. Thus, the ratio ∆F/∆H needed to obtain k (using Equation 3) can be estimated as a ratio of the corresponding slopes, i.e.,

Estimating f and h
To estimate f , we use the thermodynamic integration formula, Thus, f can be obtained using a single simulation through Alternatively, by rearranging Equation 7, we get suggesting that Ñ v κ,N * ought to differ from N * by a constant offset, − f /κ. Thus, an alternate way to estimate f from multiple simulations (with the same κ but different N * -values) would be to plot Ñ v κ,N * versus N * ; the data can then be fitted to a straight line with a unit slope, and the y-intercept used to estimate − f /κ. This procedure is functionally equivalent to averaging estimates of f obtained from different biased simulations using Equation 7.
Similarly, h can be obtained from multiple biased simulations (with the same κ) by fitting H ≡ x COM κ,N * versus N * to a straight line, and obtaining the slope. Alternatively, h can also be obtained from a single simulation using the co-variance relation where β −1 ≡ k B T , k B is the Boltzmann constant and T is the system temperature. Once both f and h have been estimated, the product of the wetting coefficient and surface tension can be readily obtained as For a surface with a known value of k, Equation 10 can then be used to obtain the vapor-liquid surface tension, γ VL . Sim-ilarly, for a liquid with a known γ VL , k can be estimated using Equation 10.

Characteristic Features of SWIPES
The setup used in SWIPES is similar in spirit to the Wilhelmy plate method, which is a commonly-used experimental technique for determining γ VL cos θ 51,52 . In the experiments, the surface of interest (the Wilhelmy plate) is immersed in the liquid, and the capillary force it experiences is measured; here we bias the fluid to wet the surface instead, and measure the corresponding free energy change per unit length ( f /h). The SWIPES approach is versatile, and can be used not only with Monte Carlo, but also with molecular dynamics simulations. It can be used to estimate both γ VL and k for fluids and surfaces that vary widely in their cohesive and adhesive interaction strengths. Moreover, SWIPES is relatively easy to implement, and because it only requires estimates of ensemble averages that converge quickly, it is also computationally efficient. Another salient feature of SWIPES is that as the number of fluid molecules in v increases, the vaporliquid interface moves along the solid surface, but no additional vapor-liquid interfacial area is created. In contrast, either vapor-liquid or vacuum-liquid interfaces are created in several computational approaches for estimating k, including the interface potential 30,31 , phantom-wall 33,34 , and dry-surface methods 35 . When carried out in the presence of periodic boundary conditions, the creation of such interfaces involves transitions between different dewetted morphologies, which can give rise to hysteresis effects, and complicate estimation of the corresponding free energies [43][44][45][46] .

C. Molecular Models and Simulation Details
To illustrate our approach for estimating γ VL and k, we study the wetting of a Lennard-Jones (LJ) solid surface by SPC/E water 50 . Following ref. 35, we study a family of surfaces composed of LJ atoms that vary in the strength of their interactions with water. The LJ atoms are placed on an FCC lattice with a spacing of 4.05 Å, and are restrained to their initial positions using harmonic springs with a spring constant of 1000 kJ/mol/nm 2 . The (111) face of the LJ surface is placed in contact with water. The LJ interaction parameters for the solid atoms are σ = 0.2629 nm and ε = 22.13 kJ/mol. The cross interaction between the solid atoms and water oxygens follows the Lorentz-Berthelot combining rule for σ , whereas the corresponding well-depth is systematically varied from ε SW = 0.001 to 2.6 kJ/mol to tune surface-water interaction strength. The cutoff distance is chosen to be 1 nm for both the LJ potential and the real space electrostatic interactions. The reciprocal space electrostatic interactions are treated using particle mesh Ewald summation with a Fourier spacing of 0.1 nm 53 . The equations of mo-tion are integrated using the leap-frog algorithm with a time step of 2 fs. The simulations are performed in the canonical ensemble, and the temperature of the system is set to 300 K using a Langevin thermostat with a coupling constant of 2 ps 54 . The SHAKE algorithm is used to constrain the internal degrees of freedom of the SPC/E water molecule 55 . All simulations were performed using GROMACS MD simulation package 56 , which was modified in-house to incorporate the biasing potentials of interest. The biasing potentials described here have also been implemented in the open-source INDUS code, which can be accessed at: http://patelgroup.seas.upenn.edu/indus.html. The INDUS code implements the INDirect Umbrella Sampling method for biasing coarse-grained particle number as an extension to the PLUMED plugin, which interfaces with several popular MD simulation packages 57 .

Biased Simulations
A solid surface consisting of 8640 particles is placed in contact with 7000 water molecules in a rectangular box with dimensions of 26.5 nm × 2.8059 nm × 14.319 nm ( Figure 1). Additionally, a repulsive wall (not shown in Figure 1) consisting of four layers of LJ particles is placed at the left edge (−x) of the simulation box to break translational symmetry and nucleate vapor in its vicinity. The wall has 1200 particles that interact with water oxygens through the LJ potential with σ = 0.29 nm and a small ε = 0.001 kJ/mol. As illustrated in Figure 1, the left edge of the observation volume is placed roughly 1 nm to the left of the solid surface of interest, and the observation volume covers the entire surface. To obtain estimates of f and h, twelve biased simulations are performed with a fixed value of κ = 0.05 kJ/mol and different N * -values that range from 4000 to 5100 in steps of 100. Each biased simulation is run for 1 ns, and the data from the first 400 ps are discarded as equilibration.
For every biased simulation,Ñ v -values are stored every 50 time steps, whereas system configurations are saved every 500 time steps, and analyzed to obtain H. System configurations are additionally used to obtain vapor-liquid interface profiles; to this end, we average the water density field over 600 configurations, and determine the iso-surface for which the averaged water density field is equal to half the bulk water density. The vapor-liquid interface profile thus obtained is fitted to a circle, and the tangent to the fitted circle at the solid surface, i.e., at the z-coordinate corresponding to the outmost layer of solid atoms, is used to obtain an estimate the contact angle, which we call θ I . Error bars are estimated by dividing the production period into 3 to 5 blocks and determining the statistical uncertainty across the block averages.

Droplet Simulations for Estimating Contact Angles
To compare the wetting coefficients obtained using our approach with those obtained from the geometry of water droplets, we also perform simulations of (infinitely long) cylindrical water droplets on the LJ surfaces described above. In these simulations, the solid surface is composed of 4200 LJ atoms, and has dimensions of 17.3607 nm × 2.8059 nm × 1.2887 nm. A box with dimensions of 17.3607 nm × 2.8059 nm × 10.0 nm is used to simulate 2000 water molecules that are placed on the surface. The system is equilibrated for 4 ns, which is followed by a production period of 6 ns. In the production stage, configurations are collected every 1000 time steps, and are analyzed to characterize the shape of the droplet. Droplet interfaces, the corresponding contact angles (θ D ), and associated error bars are obtained as described above.

D. Results and Discussions
SWIPES is used to characterize the wetting coefficients (k) for LJ surfaces with wide range of solid-water attractions (ε SW ), and to estimate the vapor-liquid surface tension of water (γ VL ). The contact angles estimated using SWIPES (θ F ) are compared with those obtained using droplet geometries (θ D ). As ε SW is increased, a transition from partial to complete wetting is observed. The computational efficiency of the method is also discussed.

Estimating wetting coefficients using SWIPES
To illustrate our method, we first estimate k for a surface with ε SW = 1.94 kJ/mol. In Figure 2a, we show the ensemble averages of the coarse-grained number of waters, Ñ v κ,N * , evaluated at different N * -values. The data are fit well by a straight line with unit slope (dashed line), confirming that f is independent of N * (Equation 8). The yintercept of the fitted line corresponds to − f /κ, enabling estimation of f = −0.143(3) kJ/mol. As shown in the Supplementary Information, estimates of f can also be obtained from individual biased simulations by using Equation 7, albeit with larger statistical uncertainties. Figure 2b shows that the position of vapor-liquid interface (H), quantified by x COM κ,N * , varies linearly with N * ; the corresponding slope is h = 8.6(2) × 10 −4 nm. As shown in the Supplementary Information, h can also be obtained from a single biased simulation by estimating the co-variance of ). However, due to the larger statistical uncertainties associated with co-variance estimation, the use of multiple biased simulations is likely to be more judicious route for estimating h. In Figure 2c, vapor-liquid interfaces identified from three biased simulations are shown with their x-coordinates shifted by x COM κ,N * ; the excellent agreement between the shifted interfacial profiles highlights the fact that biasing potential moves the water slab along the solid surface, but does not otherwise influence the shape of the interface. Thus, Figure 2c further corroborates the linear dependences of both F κ,N * and x COM κ,N * on N * . With the corresponding slopes f and h estimated, k can then be obtained from Equation 10 if γ VL is known. For the SPC/E water model used here, the surface tension has been estimated to be γ SPC/E VL = 63.6(1.5) mN/m at 300 K 58 . With this value of surface tension, we estimate the wetting coefficient for the LJ surface with ε SW = 1.94 kJ/mol to be k = 0.78 (2), and the corresponding contact angle to be 38(2) • .

Estimating vapor-liquid surface tension using SWIPES
In principle, the vapor-liquid surface tension, γ VL , of a fluid can be estimated using Equation 10 if we use SWIPES to characterize a surface (i.e., estimate f and h) for which the wetting coefficient is known. In practice, surfaces that are either fully wetting (k = 1) or non-wetting (k = −1) would work well for this purpose. However, purely repulsive surfaces are ideally suited for this purpose because they are non-wetting regardless of the fluid under consideration. Here, we use an effectively repulsive hydrophobic surface with ε SW = 0.001 kJ/mol to estimate γ VL for SPC/E water. For this surface, Ñ v κ,N * and x COM κ,N * are shown as functions of N * in Figures 3a and 3b, respectively. Using this data, we obtain estimates of f and h, which in turn enables us to estimate the vapor-liquid interfacial tension, γ VL = (1/2L) · ( f /h) = 64(2) mN/m. Our estimate is in excellent agreement with the value of 63.6 mN/m obtained in Ref. 58 using the test-area method.

LJ surfaces with a wide range of surface-water attractions
To highlight the versatility of SWIPES, we study surfaces with a wide range of solid-water interaction energies (ε SW ) that span an order of magnitude from roughly 0.2 to 2 kJ/mol. In each case, linear dependences of Ñ v κ,N * and x COM κ,N * on N * are observed; see Supplementary Information. The corresponding wetting coefficients are obtained using Equation 10, and are shown as a function of ε SW in Figure 4a. Over this range of ε SWvalues, k increases monotonically with ε SW , and varies from -0.92 to 0.91, nearly spanning the entire range of wetting coefficients possible.
As highlighted in Figure 2, the biasing potential wets the surface, but leave the curvature of the vapor-liquid interface inside v unchanged. Thus, the tangent of the vaporliquid interface at 3-phase contact can also be used to estimate the contact angle from our biased simulations. For the surface with ε SW = 1.94 kJ/mol, the contact angle is estimated to be θ I = 39(3) • (Figure 4b), and agrees well with the corresponding estimate θ F = 38(2) • obtained using  SWIPES. To compare our estimates of k with those obtained from the widely used sessile droplet method, we additionally simulate water droplets placed on the solid surfaces, and estimate the corresponding contact angles. Although spherical droplets are typically employed for this purpose, the corresponding results tend to depend on droplet size due to the influence of line tension. Following ref. 15 , we thus make use of cylindrical droplets, which are infinitely long due to our use of periodic boundary conditions. For the surface with ε SW = 1.94 kJ/mol, the droplet shape obtained from our simulations is shown in Figure 4c, and is used to estimate θ D = 40(3) • following the procedure described in section C. As shown in Figure 4a, across surfaces that span a wide range of ε SW -values, estimates of the wetting coefficients obtained from SWIPES (k F ) are found to be consistent with those obtained from interfacial profiles (k I ) and droplets (k D ). Vapor-liquid interface profiles from our biased simulations as well as geometries of cylindrical droplets for surfaces with other ε SW -values are included in the Supplementary Information.
Although the length of the 3-phase contact line remains constant in both the biased simulations and the simulations of cylindrical droplets, the corresponding contact angles (θ I and θ D ) are expected to depend on the curvature of the vapor-liquid interface, which in turn depends on the separation between the two surfaces. In contrast, our estimates of k F are not expected to be influenced by line tension or interface curvature effects, as discussed in section B. The agreement between k F and both k I and k D suggests that interface curvature are minimal (within the statistical uncertainty) for our choice of system size and inter-surface separation. As described in detail in the Supplementary Information, we studied two additional systems with smaller and larger separations, and find that k F does not depend on separation, as expected. We also estimate k I for these systems, and find any differences across different inter-surface separation are within the statistical uncertainty. Although estimates of k I and k D agree reasonably well with k F , the necessity of having to choose a contact plane introduces an inherent ambiguity in the determination of contact angles from geometric approaches [59][60][61] . Such ambiguity is expected to introduce systematic errors in contact angle estimates obtained from geometry-based methods. In contrast, the esti- Fig. 4 (a) Wetting coefficients estimated using our free energy method (k F ) are shown for LJ surfaces with a wide range of ε SW -values. These estimates are found to be consistent with the corresponding estimates obtained from the vapor-liquid interfacial profiles (k I ≡ cos θ I ), as well as those obtained from separate simulations of cylindrical droplets (k D ≡ cos θ D ). (b) For the LJ surface with ε SW = 1.94 kJ/mol, the vapor-liquid interface obtained from a biased simulation with N * = 4800 is shown (symbols). The interface is fit to a circle (red dashed line) and the contact angle θ I is obtained from the tangent to the circle at 3-phase contact. (c) The shape of a cylindrical droplet placed on an LJ surface with ε SW = 1.94 kJ/mol is shown. The droplet interface is fit to a circle (red dashed line) and the contact angle θ D is obtained from the tangent to the circle at 3-phase contact. mation of k F is free from ambiguity, and uncertainties arise only from the sampling errors in estimating Ñ v κ,N * and x COM κ,N * . Thus, highly accurate estimates of surface wetting and interfacial properties can be obtained by simply increasing the number or length of the biased simulations, making SWIPES more suitable for studying small but important changes in wettability; e.g. changes in k arising from changes in pressure, fluid composition or surface chemistry. Moreover, the determination of contact angles from droplet geometry is challenging for both the most hydrophobic and the most hydrophilic surfaces. For surfaces with contact angles greater than 140 • , the value of θ becomes increasingly sensitive to the choice of the contact plane, making its estimation from droplet geometry increasingly challenging for such hydrophobic surfaces. Conversely, for hydrophilic surfaces with contact angles below 30 • , large system sizes and long equilibration times 24 are needed for droplets to spread on the surface and achieve their equilibrium shapes.

Fully Wetting and Dewetting Surfaces
As solid-fluid interactions become more favorable, θ decreases with the system eventually transitioning from partial wetting (θ > 0 • ) to complete wetting (θ → 0 • ). Understanding the factors that govern transitions between partial and complete wetting states is of interest in diverse contexts, including the design of materials, such as polymer nanocomposites [62][63][64] . Conversely, as solid-fluid interactions are made increasingly unfavorable, the fluid is expected to dewet the solid. Although the estimation of θ from droplet geometry becomes increasingly challenging for surfaces in the vicinity of wetting and dewetting transitions, here we illustrate that SWIPES nevertheless works well in those regimes.
An LJ surface with a small value of ε SW interacts less favorably with liquid water than with its vapor phase, so that γ SV < γ SL , i.e., the surface is hydrophobic with k = (γ SV − γ SL )/γ VL < 0. As ε SW is increased, the surface interacts more favorably with both the liquid and vapor phases of water, i.e., both γ SV and γ SL decrease; however, because the denser liquid can better utilize these favorable interactions, γ SL decreases more rapidly than γ SV , and the wetting coefficient, k = (γ SV − γ SL )/γ VL increases. For a sufficiently large ε SW -value, γ SL eventually falls below γ SV , i.e., the surface becomes hydrophilic with k = (γ SV − γ SL )/γ VL > 0. Eventually, as ε SW is made large enough, any contact between solid and vapor is sufficiently unfavorable (relative to solid-liquid contact) that k → 1 and γ SV → γ SL + γ VL . In this limit, the free energetic penalty for replacing the solid-vapor interface with a solid-liquid interface and a vapor-liquid interface vanishes, and water molecules from the vapor phase adsorb onto the surface to form a film of liquid. As ε SW is increased even further, γ SL does not continue to decrease more rapidly than γ SV ; rather the two quantities continue to decrease in concert with one another, so that k plateaus at 1. The phenomenology described above is generally true regardless of the solid and the fluid wetting it 47,65 . As shown in Figure 5a, SWIPES captures such a transition from partial wetting (k < 1) to complete wetting (k = 1) as ε SW is increased, highlighting its suitability for characterizing highly hydrophilic surfaces. The ε SW -value at which the wetting coefficient, k approaches 1 is roughly ε * SW = 2.1 kJ/mol, with a further increase in ε SW leaving k unchanged. In the inset of in Figure 5a, we show a simulation snapshot of an LJ surface with ε SW = 2.4 kJ/mol (> ε * SW ) that is focused on the region near vapor-liquid-solid contact. The snapshot highlights the presence of water molecules that are adsorbed onto the solid surface, and form a liquid film which covers the surface.
In addition to the transition to complete wetting, SWIPES also captures the transition from partial wetting to dewet- ting as ε SW is decreased. For small ε SW -values, the surface is hydrophobic, and interacts more favorably with the vapor (relative to the liquid). For sufficiently small ε SW , any contact between solid and liquid is so unfavorable (relative to solid-vapor contact) that the surface dewets forming a vapor layer adjacent to the surface, in a manner that is analogous to the formation of a liquid film in the complete wetting regime. The formation of such a vapor layer results in the solid-liquid interface being replaced with a solid-vapor interface and a vapor-liquid interface, and carries a vanishing free energetic penalty as k → −1 and γ SL → γ SV + γ VL . As seen in Figure 5b, the wetting coefficient indeed approaches -1 for small ε SW . In the inset of Figure 5b, a simulation snapshot of a hydrophobic surface with ε SW = 0.05 kJ/mol is shown; the region near the surface clearly shows the presence of a vapor layer.
An important feature of SWIPES is that both the liquid and vapor phases are always explicitly present by construction. This enables water molecules to adsorb onto the surface and form a liquid film in the complete wetting regime. Similarly, the surface is readily able to dewet and form a vapor layer in the dewetting regime. Moreover, because the vapor-liquid interface is moved along the surface in SWIPES, as surface area comes in contact with liquid, it simultaneously loses contact with the vapor, enabling direct estimation of γ SV − γ SL . In contrast, methods that involve replacing the liquid adjacent to the surface with a cavity or a vacuum 33,39 instead estimate γ S,Vac − γ SL . Because γ S,Vac does not depend on ε SW , and γ SL decreases monotonically with increasing ε SW , their difference, γ S,Vac − γ SL does not plateau with increasing ε SW . Thus, such methods cannot identify the transition from partial to complete wetting unless γ S,Vac − γ SV is estimated as well 46 .

Computational Efficiency
In addition to being relative simple, SWIPES is also efficient from a computational standpoint. Here we provide a sense of the typical computational effort required to obtain accurate estimates of f and h, and thereby of k. As described in section B, the estimation of f relies on sampling ofÑ v , and an accurate estimation of its average, Ñ v κ,N * (Equation 7). Figure 6a shows the convergence ofÑ v for two biased MD simulation with N * of 4100 and 4500 for the LJ surface with ε SW = 1.94 kJ/mol. As shown in Figure 6a, the cumulative average ofÑ v − N * (black line) converges to its equilibrium value (blue dashed line) in around 200 ps 66 . Similarly, Figure 6b shows the convergence of x COM for the two biased simulations; the corresponding cumulative average converges to its equilibrium value in roughly 300 ps.
Although both f and h can also be estimated from individual biased simulations, the use of multiple simulations can yield estimates with lower statistical uncertainties. For the LJ surface with ε SW = 1.94 kJ/mol, Figure 6c shows wetting coefficient (k) estimates and uncertainties calculated using an increasing number of biased simulations; every biased simulation has a production period of 600 ps, and the N * -values of the biased simulations are separated by 100. Reasonably accurate estimates of k can be obtained with 3 to 6 biased simulations, suggesting that a total simulation time of roughly 5 ns is sufficient. This is comparable to the time needed to estimate contact angles from droplet shapes, and roughly an order of magnitude smaller than other free energy methods for estimating k, such as the procedure developed by Kanduc and Netz 47 or the dry-surface method 38 , which require 20 or more simulations, each run for 2 ns or more. We thus expect SWIPES to be comparable with or more efficient than existing methods from a computational standpoint.

E. Conclusions and Outlook
In this work, we introduce a method for the characterization of surface wetting and interfacial properties using enhanced sampling. The method, which we call SWIPES, makes use of the indirect umbrella sampling (INDUS) technique to re- versibly wet the surface of interest by biasing and systematically altering the coarse-grained number of fluid molecules in an observation volume v adjacent to the surface. SWIPES is versatile and can be used in conjunction with molecular dynamics simulations to estimate the wetting coefficient k (or equivalently the contact angle, θ ) for systems with a wide range of solid-fluid interactions. It can also be used to estimate the vapor-liquid surface tension of a fluid by estimating the work required to wet a purely repulsive solvophobic surface. Importantly, SWIPES involves displacing a vapor-liquid interface along the solid surface, but not the creation of additional vapor-liquid interfacial area. In contrast, in a number of existing free energy methods, vaporliquid interfaces spanning periodic boundaries are created; the creation of such interfaces involves transitions between different dewetted morphologies, which can give rise to hysteresis, and complicate estimation of the corresponding free energies [43][44][45] .
We use SWIPES to characterize the wetting of a family of LJ surfaces by SPC/E water. For partially wetting surfaces with a wide range of attractions for water, our predictions for the contact angles are roughly consistent with those estimated from water droplet geometries. As the attractions between the solid surface and water are increased, SWIPES is able to capture the transition from partial to complete wetting; it can also captures the dewetting transition at low attractions. Moreover, only ensemble averages ofÑ v (number of fluid molecules in v) and x COM (center of mass of fluid in the direction perpendicular to the vapor-liquid interface) are required, making SWIPES both straightforward to implement and computationally efficient; a total simulation time of roughly 5 ns is sufficient to obtain reasonably accurate estimates of k.
We note that although we have primarily discussed SWIPES in the context of surface-liquid-vapor systems, it can also be extended in a straightforward manner to other interfacial systems, such as surface-liquid-liquid or surfaceliquid-solid systems. By biasing an appropriate order parameter to control the extent to which a surface is covered by one phase relative to another, the method presented here can be generalized to characterize the preference of a surface for one liquid over another, or for a crystalline solid over its coexisting liquid 32,[67][68][69] . We further note that our approach can also be used to estimate the free energetic cost for wetting surfaces with nanoscale heterogeneities, such as chemically patterned surfaces or surfaces with nanoscale texture; for such surfaces, the notion of a contact angle itself may break down, and even if it can be defined unambiguously, estimating the contact angle from droplet geometry may be non-trivial [70][71][72][73][74][75][76] . Another important class of problems for which the thermodynamic driving force can be characterized using our approach is the infiltration of fluids into porous media 77,78 .

Supporting Information
Characterizing Surface Wetting and Interfacial Properties using Enhanced Sampling (SWIPES) The observation volume is chosen to be a cuboid. The left edge of v is placed around 1 nm to the left of the leftmost atoms in the solid surface. As shown in Figure 1 of the main text, the observation volume covers the entire solid surface, and its right edge is placed to the right of the rightmost solid atoms. Out of the two vapor-liquid interfaces that separate the liquid slab and vapor, the left interface remains well outside v (far from the surface), whereas the right interface remains well within v in all of our biased simulations. These choices enable us to ensure that as the liquid slab moves along the surface, the only physical process that occurs is the replacement of vapor-solid interfacial area by an equivalent amount of liquid-solid interfacial area, and that there are no edge effects. As the number of coarse-grained water molecules inside the observation volume (v) is increased by δÑ , the vapor-liquid interface ought to advance along the surface in x-direction by a distance δH n . However, the shape of the vapor-liquid interface is expected to be independent ofÑ ; instead, it is determined by the balance of forces at the 3-phase contact line, which is governed entirely by the corresponding interfacial tensions. Thus, δÑ and δH n ought to be related through: whereρ L (z) andρ V (z) are the coarse-grained density profiles (along the z-axis) of the liquid and vapor confined between the two solid surfaces, respectively. L is the length of the simulation box along the y-axis, and z low and z high are the lowest and highest coordinates of solid atoms that are in contact with water molecules, respectively. The corresponding slope, h n ≡ dH n /dÑ can then be expressed as: ( For a homogeneous flat solid surface, the water density profile remains unchanged regardless the position of the interface, and both z low and z high are constants. Therefore, the slope h n is a constant, i.e. H n is a linear function ofÑ . Because free energy and interface position are linearly related to one another (Equation 3 of the main text), F v (Ñ ) should also be a linear function ofÑ , making f n ≡ dF v (Ñ )/dÑ a constant. Moreover, the free energetics of the constant-Ñ ensemble, F v (Ñ ), and the free energetics of biased ensemble, F κ,N * , are related through 1,2 where U κ,N * (Ñ ) = κ 2 (Ñ − N * ) 2 is the parabolic biasing potential, and βF κ,N * v (Ñ ) = − ln P κ,N * v (Ñ ) represents the free energetics ofÑ v in the biased ensemble, with the corresponding statistics being given by P κ,N * v (Ñ ) = δ(Ñ v −Ñ ) κ,N * . Taking the derivative of Equation 3 with respect toÑ : where as discussed above, f n is a constant (independent ofÑ ). By recognizing that F κ,N * v (Ñ ) ought to have a minimum atÑ = Ñ v κ,N * , and that the correspondingÑ -derivative ought to be 0, we get: From Equation 6 of the main text, we then see that Because f n is a constant, f ≡ dF κ,N * /dN * should also be a constant. Thus, F κ,N * must be a linear function of N * . Finally, because the free energetics, F κ,N * , and interface position, H ≡ x COM κ,N * , in the biased ensemble are linearly related to one another (Equation 3 of the main text), H must also a linear function of N * , and h a constant.
For the LJ surface with SW =1.94 kJ/mol, we used umbrella sampling to obtain F v (Ñ ); i.e., we employed biased simulations with sufficient overlap in the distributions, P κ,N * v (Ñ ), between neighboring windows. We then analyzed the results using standard techniques, such as the weighted histogram analysis algorithm (WHAM) 3 to obtain both F v (Ñ ) as a function ofÑ , and F κ,N * as a function of N * . These results are shown in Figures S1a and S1b, and confirm the linear dependences of F v (Ñ ) oñ N and of F κ,N * on N * ; the corresponding slopes, f and f n also agree with one another, as expected. The simple (linear) dependences of F v (N ) and F κ,N * as well as of H n and H makes the estimation of f and h both straightforward and efficient, and lie at the heart of the simplicity and computational efficiency of SWIPES.

Estimating h from covariances
From every biased simulation, the slope, h ≡ d x COM κ,N * /dN * can be estimated using the co-variance of the fluctuations in x COM andÑ v (Equation 9 of the main text). For the LJ surface with SW = 1.94 kJ/mol, h thus obtained is shown as a function of N * in Figure S3; it is clear that the statistical uncertainties associated with the corresponding estimates of h is substantial. The average of the h-estimates obtained using the co-variance relation from the 12 biased simulations is 7.8(±1) × 10 −4 nm (dashed line in Figure S3), and agrees reasonably well with the estimate of h obtained from the slope of the line fitted to x COM κ,N * vs N * (Figure 2b of the main text). Ñ v κ,N * , are shown as a function of N * for LJ surfaces with SW -values of (a) 0.1 kJ/mol, (b) 0.5 kJ/mol, (c) 2.08 kJ/mol, and (d) 2.25 kJ/mol, respectively. Ensemble averages of water slab center of mass position along the x-axis, x COM κ,N * , are shown as a function of N * for LJ surface with SWvalues of (e) 0.1 kJ/mol, (f) 0.5 kJ/mol, (g) 2.08 kJ/mol, and (h) 2.25 kJ/mol, respectively. Dashed lines are linear fits to the simulation data (symbols). In fitting Ñ v κ,N * vs. N * , the slope is set to 1.
6 Effect of solid surfaces separation distance on the calculation of k For all the biased simulations described in the main text, the length of the simulation box in the z-direction, L z , (which sets the separation between the solid surface) was fixed to be 14.319 nm. Here we estimate k using a smaller and larger systems with L z -values of 10.023 nm and 17.182 nm, respectively. While 7000 water molecules were used in the biased simulations reported in the main text, 5000 and 12000 waters are used in the small and large systems, respectively. The solid surfaces themselves continue to be made up from 8640 atoms. In Figure S5, the ensemble averages Ñ v κ,N * and x COM κ,N * are shown as functions of N * for the small and large systems, and enable estimation of k. The contact angles thus estimated from the small and large systems using SWIPES are 40(2) • and 41(2) • , respectively, and agree with one another as well as with the system used in the main text within statistical uncertainty. Alternatively, contact angles can also be extracted from vapor-liquid interface geometries, as shown in Figure S6. The contact angles thus estimated for the small and large systems are 38(3) • and 39(3) • , respectively, and agree with one another as well as with the system used in the main text within statistical uncertainty.