Wilhelmy equation revisited: a lightweight method to measure liquid-vapor, solid-liquid and solid-vapor interfacial tensions from a single molecular dynamics simulation

We have given theoretical expressions for the forces exerted on a quasi-2D flat and smooth solid plate immersed into a liquid pool of a simple liquid. All forces given by the theory, the local forces on the top, the contact line and the bottom of the plate as well as the total force, showed an excellent agreement with the MD simulation results. We also revealed that the local forces around the bottom and top of the solid plate can be related to the SL and SV interfacial tensions \gamma_{SL} and \gamma_{SV}, and this was verified through the comparison with the SL and SV works of adhesion obtained by the thermodynamic integration. We have reached a conclusion that \gamma_{SL} and \gamma_{SV} as well as the liquid-vapor interfacial tension \gamma_{LV} can be extracted from a single equilibrium MD simulation without the computationally-demanding calculations.

We have given theoretical expressions for the forces exerted on a so-called Wilhelmy plate, which we modeled as a quasi-2D flat and smooth solid plate immersed into a liquid pool of a simple liquid. All forces given by the theory, the local forces on the top, the contact line and the bottom of the plate as well as the total force, showed an excellent agreement with the MD simulation results. The force expressions were derived by a purely mechanical approach, which is exact and ensures the force balance on the control volumes arbitrarily set in the system, and are valid as long as the solid-liquid (SL) and solid-vapor (SV) interactions can be described by meanfields. In addition, we revealed that the local forces around the bottom and top of the solid plate can be related to the SL and SV interfacial tensions γ SL and γ SV , and this was verified through the comparison with the SL and SV works of adhesion obtained by the thermodynamic integration (TI). From these results, it has been confirmed that γ SL and γ SV as well as the liquid-vapor interfacial tension γ LV can be extracted from a single equilibrium MD simulation without the computationally-demanding calculation of the local stress distributions and the TI.

I. INTRODUCTION
The behavior of the contact line (CL), where a liquid-vapor interface meets a solid surface, has long been a topic of interest in various scientific and engineering fields because it governs the wetting properties. [1][2][3][4][5] By introducing the concept of interfacial tensions and contact angle θ, Young's equation 6 is given by where γ SL , γ SV and γ LV denote solid-liquid (SL), solid-vapor (SV) and liquid-vapor (LV) interfacial tensions, respectively. The contact angle is a common measure of wettability at the macroscopic scale. Young's equation (1) was first proposed based on the wall-tangential force balance of interfacial tensions exerted on the CL in 1805 before the establishment of thermodynamics, 7 while recently it is often re-defined from a thermodynamic point of view instead of the mechanical force balance. 1 Wetting is critical especially in the nanoscale with a large surface to volume ratio, e.g., in the fabrication process of semiconductors, 8 where the length scale of the structure has reached down to several nanometers. From a microscopic point of view, Kirkwood and Buff 9 first provided the theoretical framework of surface tension based on the statistical mechanics, and molecular dynamics (MD) and Monte Carlo (MC) simulations have been carried out for the microscopic understanding of wetting through the connection with the interfacial tensions.  Most of these works on a simple flat and smooth solid surface indicated that the apparent contact angle of the meniscus or droplet obtained in the simulations corresponded well to the one predicted by Young's equation (1) using the interfacial tensions calculated through a mechanical manner and/or a thermodynamic manner, where Bakker's equation and extended one about the relation between stress distribution around LV, SL or SV interface and corresponding interfacial tension have played a key role. 21 On the other hand, on inhomogeneous or rough surfaces, the apparent contact angle did not seem to correspond well to the predicted one, 26,[35][36][37] because the pinning force exerted from the solid must be included in the wall-tangential force balance. 22 The Wilhelmy method 38 has been applied as one of the most common methods to experimentally measure the LV interfacial tension, i.e., surface tension, or the contact angle. 39 In this method, the force on a solid sample vertically immersed into a liquid pool is expressed In this work, we revisited the forces exerted on the Wilhelmy plate with non-zero thickness and derived theoretical expressions of the local forces on the CL and on the top and bottom of the plate as well as the total force on the plate. The derivations were done by a purely mechanical approach, which ensured the force balance on the arbitrarily set control volumes, and the connection to the thermodynamics was given by the extended Bakker equation. 21 We also verified the present theoretical results by MD simulations. As a major outcome of the expressions of the local forces, we will show in this article that all the interfacial tensions involved in the system, γ LV , γ SL and γ SV , can be measured from a single equilibrium MD simulation without computationally-demanding calculations.

A. MD Simulation
We employed equilibrium MD simulation systems of a quasi-2D meniscus formed on a hollow rectangular solid plate (denote by 'solid plate' hereafter) dipped into a liquid pool of a simple fluid as shown in Fig. 1. We call this system the 'Wilhelmy MD system' hereafter.
Generic particles interacting through a LJ potential were adopted as the fluid particles. The 12-6 LJ potential given by was used for the interaction between fluid particles, where r ij is the distance between the particles i at position r i and j at r j , while ǫ and σ denote the LJ energy and length parameters, respectively. This LJ interaction was truncated at a cut-off distance of r c = 3.5σ and quadratic functions were added so that the potential and interaction force smoothly vanished at r c . The constant values of c LJ 2 and c LJ 0 were given in our previous study. 19 Hereafter, fluid and solid particles are denoted by 'f' and 's', respectively and corresponding combinations are indicated by subscripts.
A rectangular solid plate in contact with the fluid was prepared by bending a honeycomb graphene sheet, where the solid particles were fixed on the coordinate with the positions of 2D-hexagonal periodic structure with an inter-particle distance r ss of 0.141 nm. The zigzag edge of the honeycomb structure was set parallel to the y-direction with locating solid particles at the edge to match the hexagonal periodicity. The right and left faces were set at x = ±x s parallel to the yz-plane, and the top and bottom faces were parallel to the xy-plane. Note that the distance between the left and right faces 2x s ≈ 1.7 nm was larger than the cutoff distance r c .
The solid-fluid (SF) interaction, which denotes SL or SV interaction, was also expressed by the LJ potential in Eq. (5), where the length parameter σ sf was given by the Lorentz mixing rule, while the energy parameter ǫ sf was changed in a parametric manner by multiplying a SF interaction coefficient η to the base value ǫ 0 sf = √ ǫ ff ǫ ss as This parameter η expressed the wettability, i.e., η and the contact angle of a hemicylindrically shaped equilibrium droplet on a homogeneous flat solid surface had a oneto-one correspondence 19,21,22 , and we set the parameter η between 0.03 and 0.15 so that the corresponding cosine of the contact angle cos θ may be from −0.9 to 0.9. The definition of the contact angle is described later in Sec. III. Note that due to the fact that the solid-solid inter-particle distance r ss shown in Table I were relatively small compared to the LJ length parameters σ ff and σ fs , the surface is considered to be very smooth, and the wall-tangential force from the solid on the fluid, which induces pinning of the CL, is negligible. 21,22 In addition to these intermolecular potentials, we set a horizontal potential wall on the bottom (floor) of the calculation cell fixed at z = z flr about 5.3 nm below the bottom of the solid plate, which interacted only with the fluid particles with a one-dimensional potential field Φ 1D flr as the function of the distance from the wall given by of T c at 90 K, which was above the triple point temperature, 45 by velocity rescaling applied to the fluid particles within 0.8 nm from the floor wall regarding the velocity components in the x-and y-directions. Note that this region was sufficiently away from the bottom of the solid plate and no direct thermostating was imposed on around the solid plate, so that this temperature control had no effects on the present results.
With this setting, a quasi-2D LJ liquid of a meniscus-shaped LV interface with the CL parallel to the y−direction was formed as an equilibrium state as exemplified in Fig. 1, where a liquid bulk with an isotropic density distribution existed above the bottom wall by choosing a proper number of fluid particles N f as shown in Fig. 2. We checked that the temperature was constant in the whole system after the equilibration run described below. Note also that in the present quasi-2D systems, effects of the CL curvature can be neglected. 14,16,19,21,22,26,46,47 The velocity Verlet method was applied for the integration of the Newtonian equation of motion with a time increment of 5 fs for all systems. The simulation parameters are summarized in Table I with the corresponding non-dimensional ones, which are normalized by the corresponding standard values based on ǫ ff , σ ff and m f .
The physical properties of each equilibrium system with various η values were calculated as the time average of 40 ns, which followed an equilibration run of more than 10 ns.

III. RESULTS AND DISCUSSION
A. Contact angle and force on the solid plate We calculated the distribution of force exerted from the fluid on the solid particles by dividing the system into equal-sized bins normal to the z-direction, where the height of the bin δz of 0.2115 nm was used considering the periodicity of the graphene structure.
We defined the average force density dξ z /dz as the time-averaged total downward (in −zdirection) force from the fluid on the solid particles in each bin divided by 2l y δz, where l y is the system width in the y-direction. Except at the top and bottom of the solid plate, dξ z /dz corresponds to the total downward force from both sides divided by the sum of surface area of both sides, i.e., the downward force per surface area. We also calculated the average SF potential energy per area u sf as well, which was obtained by substituting the downward force by the SF potential energy. This is because the time-averaged fluid density in these regions was homogeneous in the the local integral of dξ z /dz indeed gave the physical information about the force around the top and bottom parts. Note also that ξ z has the same dimension as the surface tension of force per length.
The LV interface had a uniform curvature away from the solid plate to minimize LV interface area as one of the principal properties of surface tension. Considering the symmetry of the system, the hemi-cylindrical LV interface with a uniform curvature is symmetrical between the solid plates over the periodic boundary in the x-direction. Regarding SF interface position x SF , which was different from the wall surface position x s , we defined it at the limit that the fluid could reach. With this definition, Young's equation holds for quasi-2D droplets on a smooth and flat solid surface, as shown in our previous study. 21 The x SF value was determined as x SF = 1.15 nm from the density distribution, whereas the curvature radius R was determined through the least-squares fitting of a circle on the density contour of ρ =400 kg/m 3 at the LV interface excluding the region in the adsorption layers near the solid surface. 19,21,22 We defined the apparent contact angle θ by the angle at x = x SF between the SF interface and the least-squares fit of the LV interface having a curvature χ ≡ ±1/R, with R being the curvature radius. Note that the sign ± corresponds to the downward or upward convex LV-interfaces, respectively. The relation between the SF interaction coefficient η and cosine of the contact angle cos θ is shown in Appendix A, and the following results are shown based on cos θ instead of η. Corresponding half-snapshots and density distributions are also displayed on the top. Regarding the force around the top ξ top z , it was almost zero except for cases with small contact angle. This is obvious because almost no vapor particles were adsorbed on the top of the solid plate for non-wetting cases as seen in the top panel for η = 0.03. However, in the case of large cos θ, ξ top z had non-negligible positive value, i.e., downward force comparable to ξ total z , because an adsorption layer was also formed at the SV interface as seen in the top panel for η = 0.15. In terms of the force around the contact line ξ cl z , it was positive even with negative cos θ value, meaning that the solid particle around the CL was always subject to a downward force from the fluid. On the contrary to ξ top z and ξ cl z , which were both positive, ξ bot z was negative and its magnitude increased as cos θ increased, meaning that upward force to expel the bottom side was exerted from the liquid, and that the upward force was larger for larger SL interaction η. Finally, the sum of the above three ξ total z seems to be proportional to cos θ. We will show later that it actually deviates from a simple Wilhelmy relation (3). The z-normal faces are set respectively at z = z blk V , z SV , z SL and z blk L , where z blk V and z blk L are at the vapor and liquid bulk heights, whereas z SV and z SL are set at the heights of SV and SL interfaces, respectively as shown in Fig. 3 at which dξ z /dz = 0 is satisfied. These heights can be set rather arbitrary as long as the above conditions are satisfied. We define the forces from solid to liquid by F top z , F cl z and F bot z on the top, middle and bottom CVs, respectively.
In addition, we also categorize the right-half of the solid plate into top, middle and bottom parts shown with dark-yellow, blue, and red solid lines, respectively with z SV and z SL as the boundaries as shown in Fig. 4. where forces ξ top z , ξ cl z and ξ bot z in the z-direction are exerted from the fluid, respectively. Specifically note that ξ cl z also includes the forces from the top and bottom parts of the solid, whereas ξ cl z includes the forces from the top and bottom CVs. In other words, the force between the middle solid part and middle fluid CV is in action-reaction relation, but F cl z and ξ cl z include different extra forces above. This will be described more in detail in the following.

Capillary force ξ cl z around the contact line based on a mean-field approach
We start from formulating the wall tangential force on the solid particles ξ cl z on the right face of the solid plate. Taking into account that the solid is supposed to be smooth for the fluid particles because the interparticle distance parameters σ ff and σ sf are sufficiently large compared to r ss between solid particles, ξ cl z can be analytically modeled by assuming the mean fields of the fluid and solid. The mean number density per volume of the fluid is given as a function of the two-dimensional position (z f , x f ) of the fluid, whereas a constant mean number density per area ρ s A of the solid is used considering the present system with a solid plate of zero-thickness without volume; however, the following derivation can easily be extended for a system with a solid with a volume and density per volume in the range x ≤ x s as long as the density is independent of z s . We start from the potential energy on a solid particle at position (x s , y s , z s ) due to a fluid particle at (x f , y f , z f ) given by Eq. (5). We define in the following. Assuming that the fluid particles are homogeneously distributed in the y-direction with a number density ρ f V (z f , x f ) per volume, the mean potential field from an infinitesimal volume segment of dz f × dx f on the solid particle is defined by using This schematic is shown in the inset of Fig. 5. Then, the local tangential force f s z (z ′ f , x ′ f )dz f dx f dz s exerted on an infinitesimal solid area-segment of dz s from the present fluid volume-segment is given by: where denotes the tangential force density on the solid. Note that dx f and dx ′ f are identical because x s is a constant.
Since Φ LJ (r) is truncated at the cutoff distance r c in the present case, for ) as a function of x ′ f denotes the cutoff with respect to z ′ f . Indeed this cutoff is not critical as long as φ (z ′ f , x ′ f ) quickly vanishes with the increase of r, but we continue the derivation including the cutoff for simplicity. With the definition of x SF as the limit that the fluid could reach, it follows that In addition, considering that φ(z ′ f , x ′ f ) is an even function with respect to z ′ f , i.e., it follows for the mean local potential φ that and where Eq. (9) is applied for the latter, which corresponds to the action-reaction relation between solid and fluid particles under a simple two-body interaction, i.e., holds for the tangential force density on the fluid f f z . Based on these properties, we now derive the analytical expression of ξ cl z as the triple integral of the local tangential force f s z in Eq. (12) around the CL, where the fluid density ρ f V decreases with the increase of z f within a certain range. Let this range be and let ρ f V outside this range be given as a unique function of x f by as shown in Fig. 5. Then, ξ cl z is expressed by as the triple integral of the force density f s z in Eq. (13), where the integration range of the double integral regarding z f and z s corresponds to the region filled with blue in Fig. 5.
To obtain the double integral as the square brackets in Eq. (22) for the blue-filled region in Fig. 5, we calculate at first that in the region surrounded by the solid-blue line, add those in the vertically-hatched regions, and subtract those in the horizontally-hatched regions.
by using Eq. (16). Indeed, from Eq. (19), the reaction force −F cl z from solid on the fluid around the CL in the blue-dotted line in Fig. 4 is obtained by further integrating Eq. (23) with respect to x f , i.e., The final equality means that no tangential force acts on the fluid there as mentioned in our previous study. 21 Regarding the bottom-left vertically-hatched region in Fig. 5, the double integral is where φ(z c , x ′ f ) = 0 and Eq. (17) is used for the 2nd equality. This region physically corresponds to the interaction between blue solid part and fluid in the red-dotted part in Fig. 4. For the bottom-left horizontally-hatched region in Fig. 5, it follows that This region corresponds to the interaction between red solid part and fluid in the blue-dotted part in Fig. 4. Hence, the net force due to the double integral in the bottom-left hatched regions in Eqs. (25) and (26) with also integrating in the x f -direction, which we define by As a physical meaning, u SL represents the SL potential energy density i.e., potential energy per SL-interfacial area at the SL interface away from the CL and from the bottom of the solid plate.
Regarding the top-right hatched regions, the net force results in −u SV with the SV potential energy density area given by which can be derived in a similar manner. Thus, it follows for the force −ξ cl z from the fluid on the solid around the CL that therefore, by using F cl z = 0 in Eq. (24), is derived as the analytical expression of ξ cl z , where the final expression is appended considering that the potential energy densities u SL and u SV are both negative. Figure 6 shows the dependence of the SL and SV potential energy density u SL and u SV , respectively as the potential energies per interfacial area, on the cosine of the contact angle cos θ, and comparison between the force on the solid around the CL ξ cl z and difference of potential energy density −u SL + u SV . Very good agreement between ξ cl z and −u SL + u SV is observed within the whole range of the contact angle, and this indicates that Eq. 30 is applicable for the present system with a flat and smooth surface. It is also qualitatively apparent from Eq. (30) that ξ cl z is positive regardless of the contact angle because the SF potential energy is smaller at the SL interface than at the SV interface. It is also interesting to note that for the very wettable case with large cos θ, i.e., large wettability parameter η, ξ cl z decreased with the increase of cos θ. This can be explained as follows: the change can also be included in the inter-molecular force contribution, but only the FF interaction as the internal force is taken into account as the stress, and SF contribution is considered as an external force in this study. 11,21,22,51,52 With this setting, the stress is zero at the SF boundary for all CVs because no fluid particle exists beyond the boundary to contribute to the stress component as the kinetic nor at inter-molecular interaction contribution. Hence, the force balance on each CV containing only fluid is satisfied with the sum of the stress surface integral and external force from the solid. The force balance on the red-dotted CV in Fig. 4 in the z-direction is expressed by with the stress contributions from the bottom and top and external force in the RHS, respectively, by taking into account that τ xz = 0 on the x-normal faces at x = 0 and x = x end due to the symmetry, and also that the stress at the SF interface is zero. Similarly, the force balance on the blue-dotted CV and dark-yellow-dotted CV in Fig. 4 in the z-direction are expressed by and − x end respectively.
By taking the sum of Eqs. (34), (35) and (36), and inserting Eq. (33), it follows for ξ total z that ξ total Since the bottom face of the red-dotted CV and top face of the dark-yellow-dotted CV in Fig. 4 are respectively set in the liquid and vapor bulk regions under an isotropic static pressure p blk L , and p blk V given by and the 1st and 2nd terms in the RHS of Eq. (37) write and Thus, Eq. (37) results in a simple analytical expression of Furthermore, by applying the geometric relation with χ being the LV interface curvature and the Young-Laplace equation for the pressure difference in Eq. (42): which hold irrespective of whether the LV-interface is convex downward or upward, it follows for Eq. (42) as another analytical expression of ξ total z that ξ total   42) is simply the force balance to be satisfied for equilibrium systems. Regarding the pressure, p blk V is almost constant, which corresponds to the saturated vapor pressure at this temperature. In addition, a linear relation between p blk L −p blk V and cos θ can be observed, and this indicates that the Young-Laplace equation (44) is applicable in the present scale. We evaluated γ LV from this relation with the least-squares fitting, and the resulting value was γ LV = 9.79 ± 0.23 × 10 −3 N/m with x SF = 1.15 nm and x end = 7.5 nm, which was indeed close to the value obtained by a standard mechanical process. 18 The standard Wilhelmy equation (3) using this value is also shown in Fig. 7, indicating that γ LV would be overestimated with this standard Wilhelmy equation (3) in a small measurement system like the present one.
Finally, we derive the analytical expression of the local force ξ bot z and ξ top z . For the derivation of ξ bot z , we apply the extended Bakker's equation for the SL relative interfacial tension 21,22 for the 2nd term in the LHS of Eq. (34), where γ SL − γ S0 is the SL interfacial tension relative to the interfacial tension between solid and fluid with only repulsive interaction (denoted by "0" to express the solid surface without adsorbed fluid particles). Then, it follows that By inserting Eqs. (31), (40) and (47) into Eq. (34), the analytical expression of ξ bot Similary, by applying the Extended Bakker's equation for the SV interfacial tension 21,22 to Eq. (36) with Eq. (32), the analytical expression of ξ top To verify Eqs. (48) and (50), we compared the present results with ξ bot z and ξ top z calculated using the corresponding SL and SV works of adhesion W SL and W SV obtained by the thermodynamics integration (TI) with the dry-surface scheme. 21,28 The calculation detail is shown in Appendix B. By definition, the SL and SV interfacial tensions γ SL and γ SV are related to W SL and W SV by and respectively, where the approximation γ L0 ≈ γ LV for the interfacial tension γ L0 between liquid and vacuum is used in Eq. (51), and γ V0 is set zero in the final approximation in Eq. (52). Note that γ L0 or γ LV is included in W SL . From Eqs. (51) and (48), and from Eqs. (52) and (50), ξ bot z and ξ top z are respectively rewritten by and The error bar for ξ bot z using W SL in blue comes from the evaluation of γ LV from p blk L and p blk V in Fig. 7.

C. Discussion
We list the key issues for the further application of the present expression in the following. First, Eqs. (34), (35) and (36) are about the force balance and should be satisfied in equilibrium systems without any restrictions. In addition, Eqs. (29), (31) and (32) are about the relation between the solid-fluid and fluid-solid forces and should hold as long as the solid plate can be decomposed into the three parts without the interface overlapping. At both SL and SV interfaces, which are between the CL and the plate bottom and between CL and the plate top respectively, a quasi-one-dimensional density distribution with ∂ρ/∂z = 0 can be assumed and one can apply the mean-field approach described in Sec. III B 2. Furthermore, Eqs. (46) and (49) Considering the above discussion, we summarize the procedure to extract the wetting properties. In a single Wilhelmy MD simulation, we can calculate 1. Force ξ top z , ξ cl z and ξ bot z on three parts of the solid from the force-density distribution dξ z /dz in the surface-tangential direction, 2. SF potential energy densities u SL and u SV on solid per area at SL and SV interfaces, respectively from the distribution of the potential energy density u sf , 3. Bulk pressures p blk V and p blk L measured on the top and bottom of the system, and 4. Contact angle θ from the density distribution.
From these quantities the following physical properties can be obtained: a. SL relative interfacial tension γ SL − γ S0 from ξ bot z , u SL , x SF and p blk L using Eq. (48), b. SV relative interfacial tension γ SV − γ S0 from ξ top z , u SV , x SF and p blk V using Eq. (50), c. LV interfacial tension γ LV from p blk V , p blk L , x SF , the system size x end and the contact angle θ using Eq. (44) , and d. Pinning force F cl z from Eq. (29) to be added to Young's equation, which is zero in the case of flat and smooth solid surface.
Related to the above procedure, it should also be noted that, surprisingly, the microscopic structure of the bottom face does not have a direct effect on the force ξ bot z . This is similar to buoyancy given by the 3rd term of the RHS of Eq. (2), which depends on the volume V immersed into the liquid and is not directly related to the microscopic structure.
Finally, we compare the present analytical expression of the contact line force ξ cl z with an existing model by Das et al. 41 , which states This model is derived based on the assumption that the densities of the liquid and vapor are constant at bulk values even close to the solid interface: the so-called sharp-kink approximation. This is similar to the interface of two different solids whose densities and structures do not change upon contact. Even under this assumption, the force ξ cl z on solid around the CL is expressed by Eq. (30) as the difference between the SL and SV potential energy densities u SL and u SV as well. 41 The difference arises for the works of adhesion. Under the sharpkink approximation, it is clear that the works of adhesion required to quasi-statically strip the liquid and vapor off the solid surface are equal to the difference of solid-fluid potential energies after and before the procedure, i.e., which indeed results in Eq. (56) with Eqs (51) and (52). However, the density around the solid surface is not constant as shown in the density distribution in Fig. 2, and the difference of W SL and W SV is not directly related to the SL and SV potential energy densities u SL and u SV as in Eq. (57). In other words, the fluid can freely deform and can have inhomogeneous density in a field formed by the solid at the interface to minimize its free energy at equilibrium, and this includes the entropy effect in addition to u SL and u SV as parts of the internal energies. 33

IV. CONCLUSION
We have given theoretical expressions for the forces exerted on a Wilhelmy plate, which we modeled as a quasi-2D flat and smooth solid plate immersed into a liquid pool of a simple liquid. By a purely mechanical approach, we have derived the expressions for the local forces on the top, the contact line (CL) and the bottom of the plate as well as the total force on the plate. All forces given by the theory showed an excellent agreement with the MD simulation results.
In particular, we have shown that the local force on the CL is written as the difference of the potential energy densities between the SL and SV interfaces away from the CL but not generally as the difference between the SL and SV works of adhesion. On the other hand, we have revealed that the local forces on the top and bottom of the plate can be related to the SV and SL works of adhesion, respectively. As the summation of these local forces, we have obtained the modified form of the Wilhelmy equation, which was consistent with the overall force balance on the system. The modified Wilhelmy equation includes the cofactor taking into account the plate thickness, whose effect can be significant in small systems like the present one.
Finally, we have shown that with these expressions of the forces all the interfacial tensions γ SL and γ SV as well as γ LV can be extracted from a single equilibrium MD simulation without the computationally demanding calculation of the local stress distributions and the thermodynamic integrations. in the adsorption layers near the solid surface. Figure 9 shows the relation between the SL interaction parameter η and the apparent contact angle θ. The contact angle cosine cos θ monotonically increased with the increase of η, and a unique relation can be obtained between the two for the present range of η.

Appendix B: Thermodynamic integration (TI) with the dry-surface scheme
We calculated the solid-liquid (SL) and solid-vapor (SV) works of adhesion W SL and W SV , respectively, by the thermodynamic integration (TI) 53 through the dry-surface (DS) scheme 28 to compare with the relative SL and SV interfacial tensions obtained in the present Wilhelmy MD systems. Details of the DS scheme were basically the same as in our previous study. 21 In the systems shown in Fig. 10, the liquid or vapor was quasi-statically stripped off from the solid surface fixed on the bottom of the coordinate system, which had the same periodic honeycomb structure as the solid plate in the Wilhelmy MD system. The work of adhesion was calculated as the free energy difference after and before the above procedure, where the coupling parameter for the TI was embedded in the SF interaction parameter in the DS scheme.
For the calculation of W SL , a SL interface was formed between the liquid and bottom solid as shown in Fig. 10 (a) with wettability parameter η corresponding to the Wilhelmy MD system. Periodic boundary condition was employed in the x-and y-directions tangential to the solid surface. In addition, we set a piston at z = z pis above the liquid to attain a constant pressure system. By allocating sufficient number of fluid particles N f and by setting the pressure p set above the vapor pressure, a liquid bulk with a constant density was formed between the solid wall and piston. We used 3000 fluid particles, and the system size was set as shown in Fig. 10 (a). We also controlled the temperature of the fluid particles within 0.8 nm from the top piston regarding the velocity components in the x-and y-directions at We embedded a coupling parameter λ into the SF interaction potential given in Eq. (5) as and we obtained multiple equilibrium systems with various λ values with 0 ≤ λ < 1 to numerically calculate the TI described below. Each system was obtained after a preliminary equilibration of 10 ns, and the time average of 20 ns was used for the analysis.
The work of adhesion W SL is defined by the minimum work needed to strip the liquid from the solid surface per area under constant NpT , and it can be calculated by the TI along a reversible path between the initial and final states of the process. In the present DS scheme, this was achieved by at first forming a SL interface, and then by weakening the SF interaction potential through the coupling parameter. We obtained equilibrium SL interfaces with discrete coupling parameter λ varied from 0 to 0.999. Note that the maximum value of λ was set slightly below 1 to keep the SF interaction to be effectively only repulsive.
This value is denoted by 1 − hereafter. The difference of the SL interfacial Gibbs free energy ∆G SL ≡ G SL | λ=1 − − G SL | λ=0 between systems at λ = 0 and λ = 1 − under constant NpT was related to the difference in the surface interfacial energies as where the vacuum phase was denoted by subscript '0' and γ S0 and γ L0 were the solid-vacuum and liquid-vacuum interfacial energies per unit area. Note that γ L0 was substituted by the liquid-vapor interfacial tension γ LV in the final approximation considering that the vapor density was negligibly small. Using the NpT canonical ensemble, the difference of the SL interfacial Gibbs free energy ∆G SL in Eq. (B2) was calculated through the following TI: where H was the Hamiltonian, i.e., the internal energy of the system and N w was the numbers of wall molecules. The ensemble average was substituted by the time average in the simulation, and was denoted by the angle brackets. Note that to obtain ∆G SL , the work exerted on the piston Ap set ( z p | λ=1 − − z p | λ=0 ) was subtracted from the change of the Gibbs free energy of the system ∆G including the piston in Eq. (B4).
For the calculation of the SV work of adhesion W SV , we investigated the interfacial energy between saturated vapor and corresponding solid surface set on the bottom of the simulation cell by placing an additional particle bath on the top as shown in Fig. 10 (b). The setup regarding the periodic boundary conditions employed in x-and y-directions, temperature control and placement conditions for the solid surface were the same as the SL system, whereas the particle bath was kept in place by a potential field at a fixed height sufficiently far from the solid surface. This potential field mimicked a completely wettable surface with an equilibrium contact angle of zero with the present potential parameters, i.e., a liquid film was formed on the particle bath. With this setting, a solid-vapor interface with the same density distribution as that in the Wilhelmy MD system was achieved. We formed multiple equilibrium systems with various values of the coupling parameter λ with the same recipe as the SL systems.
Similar to the calculation of W SL , the SV interface at λ = 0 was divided into S0 and V0 interfaces at λ = 1 − as shown in Fig. 10 (b), while the calculation systems for W SV were under constant NV T . Thus, the solid-vapor work of adhesion W SV was given by the difference of the Helmholtz free energy ∆F per unit area, and was related to the difference in the surface interfacial energy as where γ V0 was set zero in the final approximation. Using the NV T canonical ensemble, ∆F in Eq. (B5) was calculated through the TI as: (B6) Figure 11 shows the SL and SV works of adhesion W SL and W SV calculated by the TI as a function of the solid-fluid interaction coefficient η. These values were used for the results shown in Fig. 8 through η-cos θ relation in Fig. 9.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.