On the equilibrium contact angle of sessile liquid drops from molecular dynamics simulations

We present a new methodology to estimate the contact angles of sessile drops from molecular simulations by using the Gaussian convolution method of Willard and Chandler [J. Phys. Chem. B 114, 1954-1958 (2010)] to calculate the coarse-grained density from atomic coordinates. The iso-density contour with average coarse-grained density value equal to half of the bulk liquid density is identified as the average liquid-vapor (LV) interface. Angles between the unit normal vectors to the average LV interface and unit normal vector to the solid surface, as a function of the distance normal to the solid surface, are calculated. The cosines of these angles are extrapolated to the three-phase contact line to estimate the sessile drop contact angle. The proposed methodology, which is relatively easy to implement, is systematically applied to three systems: (i) a Lennard-Jones (LJ) drop on a featureless LJ 9-3 surface; (ii) an SPC/E water drop on a featureless LJ 9-3 surface; and (iii) an SPC/E water drop on a graphite surface. The sessile drop contact angles estimated with our methodology for the first two systems are shown to be in good agreement with the angles predicted from Young's equation. The interfacial tensions required for this equation are computed by employing the test-area perturbation method for the corresponding planar interfaces. Our findings suggest that the widely adopted spherical-cap approximation should be used with caution, as it could take a long time for a sessile drop to relax to a spherical shape, of the order of 100 ns, especially for water molecules initiated in a lattice configuration on a solid surface. But even though a water drop can take a long time to reach the spherical shape, we find that the contact angle is well established much faster and the drop evolves toward the spherical shape following a constant-contact-angle relaxation dynamics. Making use of this observation, our methodology allows a good estimation of the sessile drop contact angle values even for moderate system sizes (with, e.g., 4000 molecules), without the need for long simulation times to reach the spherical shape.


I. INTRODUCTION
The contact angle, the angle between a two-fluid interface and a solid substrate, is crucial in the characterization and quantification of the wetting properties and wettability characteristics of fluid-solid systems. [1][2][3][4] It is the result of the balance between the interfacial tensions of the different phases involved and is expressed by Young's equation. 5 When one of the two fluids is liquid and the other is its own vapor, Young's equation is written as where θ Y is Young's contact angle and γ sv , γ sl , and γ lv are the solid-vapor, solid-liquid, and liquid-vapor (LV) interfacial tensions, respectively. The contact angle then plays the role of a boundary condition via Young's equation. It is a) Electronic mail: s.kalliadasis@imperial.ac.uk b) Electronic mail: a.galindo@imperial.ac.uk noteworthy that the equation reflects the mechanical equilibrium in a direction parallel to the solid surface at the threephase contact line (the line formed by the intersection between the two-fluid interface and the solid surface). But there is also a relation in the normal direction which, in the framework of the statistical mechanics of classical fluids, can be shown to be closely connected to the concept of disjoining pressure (e.g., Refs. 4 and 6-8).
Two popular approaches to obtain the contact angles from molecular dynamics (MD) simulations have been followed over the years. In one of them, the contact angle is calculated by applying Young's relation [Eq. (1)] with the interfacial tension values obtained by simulating the corresponding planar interfaces. Early studies followed the mechanical route 9 involving the calculation of the components of the pressure tensor to compute the required interfacial tension values (which in turn are used in the calculation of the contact angle from Young's relation). In addition to the mechanical route, there have also been a number of new techniques proposed to calculate interfacial tensions, such as the test-area perturbation method based on the thermodynamic definition of interfacial tension [10][11][12] and free-energy-based methods which rely on transition matrix Monte Carlo simulations. [13][14][15] In the phantom-wall method, 16,17 the work of adhesion (W = γ lv + γ sv γ sl ) is computed using thermodynamic integration and the contact angle is calculated in conjunction with the LV interfacial tension.
Although this first approach has been widely used for fluids on planar substrates, its application for substrates with topographical and chemical heterogeneities is non-trivial. Moreover, to study three-phase contact line features such as line tension, dynamic wetting, and spreading effects, simulation of planar interfaces is not helpful. For these purposes, a second approach consisting of a direct simulation involving stabilisation of a sessile drop of a given fluid on the substrate of interest, and analysing the interface to extract the contact angle, is the preferred route. Once the droplet is stabilised, the local fluid density is calculated by slicing the simulation box into bins, counting the number of particles in each bin, and calculating the ratio of the number of particles in a given bin to the volume of the bin. The collection of interfacial bins, bins with average density equal to half of the bulk liquid density (iso-density contour), constitutes the LV interface. In general, the droplet radius [r(z)] as a function of normal distance from the solid surface (z) is then fitted to a circle, [18][19][20][21] in effect assuming the LV interface to be a spherical cap. The contact angle is calculated as the angle of the tangent of the fitted circle with the solid surface at the three-phase contact line. In this approach, it is a common practice to ignore the simulation data below a pre-defined cut-off distance to avoid the influence of strong fluid density fluctuations close to the substrate on the contact angle.
Not surprisingly, there have been a large number of studies considering how the two approaches outlined above complement each other, for instance, in the case of Lennard-Jones (LJ) sessile drops. Saville 22 in his pioneering study of 1977 computed the contact angles for a LJ fluid on a LJ 9-3 planar surface using MD simulations and reported disagreement between the values obtained from the two methods leading to the conclusion that Young's equation [Eq. (1)] was not valid for nanodrops. It is important to note that Saville visually approximated the tangent to the iso-density contour close to the wall for the contact angle and that the maximum system size considered in his study was limited to 1205 particles. Interestingly, in a later study, Nijmeijer et al. [23][24][25] repeated the comparison and concluded that Young's equation was valid in the case of nanodrops. But instead of the visual approximation of the tangent taken by Saville, Nijmeijer et al. fitted the LV interface to a circle and evaluated the tangent from the equation of the circle, thus obtaining the contact angle. In essence, the difference in the results of the two studies can be traced back to the level of rigor in the estimation of the tangent to the interface. In recent years, independent studies by Ingebrigtsen and Toxvaerd, 26 Weijs et al., 20 and Peng, Birkett, and Nguyen 21 have confirmed the validity of Young's equation [Eq. (1)] using calculation of contact angles of nanodrops via MD simulations. Furthermore, in addition to validating Young's equation, the agreement between the two approaches also serves as a benchmark of the various methods used for sessile drops and concomitant contact angle calculation.
Despite the considerable attention that nanodrop simulations have received for several decades and even though analysis of LJ sessile drops using the spherical-cap approximation is now well established, there are still a degree of uncertainty and conflicting reports in its application for liquids with anisotropic interactions, specifically in the case of water drops. Hautman and Klein 27 in their study of microscopic wetting phenomena on different substrates fitted the LV interface of the water sessile drop to a sphere: the radius of the droplet [r(z)] vs the normal distance from the solid surface (z) was fitted to a circle-which is what we have referred to as the spherical-cap approximation. Later studies involving water sessile drops 28-33 adopted the methodology of Hautman and Klein to estimate sessile drop contact angle values. In these studies, a variation of the contact angle values for different drop sizes was noted and was attributed to the effect of line tension. The macroscopic contact angle, i.e., the contact angle of an infinite size droplet, was estimated by extrapolating the contact angle values of droplets of different sizes. The effect of line tension has also been reported for LJ drops 20,21 but is not as significant as that for water drops. 29,31,34 Ingebrigtsen and Toxvaerd 26 cautioned that line tension might be inclusive of the other effects, such as curvature dependence of the LV interfacial tension, but also of shortcomings of the particular ways the contact angle is estimated from the drops. In fact, the subject of line tension and its effects on contact angle remain to this day an active area of debate and research.
Setting aside the possibility of line tension effects, let us also briefly review studies centered on water which have considered other approaches to estimating the sessile drop contact angles. Giovambattista, Debenedetti, and Rossky 35 and Li and Zeng 36 relaxed the assumption of sphericity, adopting instead a generic quadratic fit of the r(z) vs z profile, while Zhang et al. 37 used a parabola. Unfortunately, the reasons for the choice of non-circular fits to the r(z) vs z data and a comparison to the corresponding macroscopic (Young's) value are not provided in these studies. More recently, Santiso, Herdes, and Müller 34 have followed a different approach to extract the water sessile drop contact angle, with emphasis on addressing the large fluctuations observed in the shape of the LV interface. These authors explored large drops, of ∼50 nm diameter (∼9 × 10 5 beads), using a coarse-grained water model based on a Mie potential, 38 in which a single bead represents two water molecules. They used surface meshing techniques to define a geometric surface by treating the interfacial fluid molecules as a point cloud data set ("local surface") and calculated the contact angle as the arc-cosine of the component of the averaged normal to the interface along the normal to the solid surface. The authors concluded that one may need to simulate sessile drops with approximately half a million water molecules to recover the macroscopic contact angle. Along the same lines, Khalkhali et al. 39 applied the Quickhull 40 approach to triangulate the sessile drop LV interface and obtained a probability distribution of the contact angle values from the normals to the triangles, but with a significant spread in the contact angle values. In addition, Shih et al. 41 and Włoch et al. 42 have reported simulation times of 100 and 50 ns, respectively, to be required for the contact angle to reach equilibrium for a system of 2000 water molecules on a graphite surface, which are very different to the time scales usually encountered in previous literature (5-10 ns).
Here, we propose a new methodology to estimate sessile drop contact angles and validate it by comparison with Young's equation. The interfacial tensions required in the equation are computed by employing the test-area perturbation method [10][11][12] for the corresponding planar interfaces. The angles between the unit normal vectors to the average LV interface and unit normal vector to the solid surface, as a function of the distance normal to the solid surface, are calculated and extrapolated to the three-phase contact line. Not only does the methodology not require the presumption of the spherical-cap approximation of the LV interface but it also allows us to scrutinize the geometry of the interface by looking at the cosines of the angles as a function of the normal distance from the solid surface (with a linear trend reflecting the spherical shape of the LV interface). (For comparison purposes, the sphericalcap approximation of the LV interface is adopted for LJ drops.) The density calculation is based on the Gaussian convolution technique proposed by Willard and Chandler, 43 which makes it easier to compute the normal at a point on the LV interface; details are discussed in Sec. III A. The spherical-cap approximation is outlined in Sec. III B. One of our objectives is to shed light on the ambiguity regarding system sizes and simulation times of drops, in particular, water drops. For this purpose, we consider three prototypical systems: (i) a LJ drop on a LJ 9-3 surface, (ii) a water drop on a LJ 9-3 surface, and (iii) a water drop on a graphite surface. The extended single point charge (SPC/E) 44 model is used for water with the forcefield details given in Sec. II B. Different contact angle scenarios are considered by varying the fluid-substrate interaction strength.
Our study reveals that water sessile nanodrops can take a long simulation time to equilibrate to a spherical shape, of the order of 100 ns starting from a lattice configuration. As may be expected, this time scale varies with different initial configurations as well as different fluid-substrate interactions, and hence it is not easy to a priori estimate the simulation time required for the equilibration. In this context, the sphericalcap approximation has to be used with care. With sufficiently long simulation time, so that the water sessile drop reaches the spherical shape, the spherical-cap approximation yields good estimates of the contact angle values, without the need to resort to large system sizes (e.g., as large as half a million molecules). We also find that the contact angle relaxes much faster than the typical time scale required for a sessile drop to adopt the spherical shape (i.e., if the spherical-cap approximation is to be used, a longer simulation time is usually needed starting from a lattice configuration). Moreover, it is shown that sessile drops follow a constant-contact-angle relaxation dynamics toward the spherical shape.

II. MOLECULAR MODEL AND SIMULATION DETAILS
We consider two fluids, Lennard-Jones (LJ) fluid and SPC/E water, and two surfaces, a featureless LJ 9-3 surface (for both LJ fluid and water), and a graphite surface for water. Further details are provided in the following.

A. Lennard-Jones fluid
The truncated LJ potential 45 is given by , and the fluid-LJ 9-3 surface potential (function only of the distance normal to the solid surface) is given by where reduced units defined in terms of the fluid-fluid collision diameter σ ff and fluid-fluid potential well depth ε ff are used. Here r * = r/σ ff is the distance between two fluid molecules, ε * sf = ε sf /ε ff is the solid-fluid potential well depth, is the fluid density, T * = Tk B /ε ff is the temperature, and k B is the Boltzmann constant. The reduced fluid-fluid potential well depth is ε * ff = 1 and the fluid-fluid LJ interaction and solidfluid LJ interaction are cut-off at r * c = z * c = 5, without any tail correction.

B. Water
In our simulations, water is represented using the rigid SPC/E model. 44 This is a three site model with a LJ site centered on the oxygen (O) atom, with the potential characterized by well depth ε OO = 0.155 kcal/mol and collision diameter σ OO = 3.166 Å. In this model, the oxygen atom is associated with a charge q O = 0.8476 e (1 e = 1.602 × 10 19 C) and the two hydrogen atoms are associated with q H = +0.4238 e, placed at a distance r OH = 1.000 Å from the oxygen atom forming an angle HOH = 109.47 • . The total interaction potential is a combination of a truncated-LJ potential between oxygen atoms in the water molecules and the electrostatic potential, i.e., where q i and q j are the charges on atoms i and j, respectively, and 0 = 8.854 × 10 12 F/m is the permittivity of the vacuum. The LJ interaction u LJ (r) between two oxygen atoms, separated by distance r, is given by We do note however, that despite the academic value of SPC/E as a prototypical molecular model for water, if direct comparisons with experimental data are envisaged, other options are advisable, such as by Vega and de Miguel. 46

C. Simulation details
In the simulations involving water, the fluid-fluid and the fluid-solid LJ interactions are cut-off at a distance of 12.5 Å, without the tail correction. The real-space electrostatic interaction is also cut-off at 12.5 Å, and the long-range electrostatic interactions are handled using the particle-particle-particle mesh (pppm) method. 47 The SHAKE algorithm 48 is used to maintain the rigidity of the water molecules. In the case of water on graphite, two graphene layers separated by 3.420 Å in the z-direction (periodic in the xand the y-direction) are prepared as a substrate. As in previous studies, 28,29,31 the carbon atoms in the graphite are frozen. The interaction between the oxygen atom in water and a carbon atom in graphene is modeled via a 126 LJ potential characterised by ε CO = 0.094 kcal/mol and σ CO = 3.190 Å, and the same cutoff of 12.5 Å is applied, without any tail correction.
All simulations are carried out in the canonical ensemble (constant number of particles, volume, and temperature) using the Nosé-Hoover thermostat in LAMMPS, an open source MD software. 49 In bulk fluid simulations, periodic boundary conditions are applied in all three directions. Simulations in the presence of a substrate are carried out with periodic boundary conditions only in the xand the y-direction and nonperiodicity in the z-direction. Specifically, at the top side of the non-periodic dimension, we place a reflective wall (velocity of a fluid molecule is simply reversed) for the case of a LJ fluid or a repulsive wall (a LJ 9-3 surface with the fluid-substrate potential cut-off distance set as σ oo ) for the case of water. Electrostatics in slab geometries (periodic in the xand the y-direction and non-periodic in the z-direction) are handled using the method proposed by Yeh and Berkowitz 50 and implemented in LAMMPS. Different contact angle scenarios are generated through different values of ε sf and keeping σ sf equal to σ ff . For water on graphite, all long-range electrostatic interactions are pre-defined from bulk studies, 44 while the interaction parameters are taken from Werder et al., 28 which were optimized to reproduce the macroscopic contact angle of water on graphite, 86 • . Although Werder et al. truncated electrostatic interactions at 10 Å, we considered long-range electrostatic interactions in our study. Let us also note here the recent experimental findings concerning the effects of defects on the intrinsic wettability of graphite 51 and also rigorous forcefield developments through parameterizing water-graphene nonbonded interactions from ab initio calculation data 52 or using the work of adhesion. 53 For all sessile drop simulations, the initial configuration consists of fluid molecules arranged in a lattice on a solid surface. In the case of a LJ sessile drop on a LJ 9-3 surface, we consider 5072 particles in a simulation box with dimensions 60 × 60 × 50. The system is equilibrated for 5 × 10 5 time steps until a sessile drop is formed and is followed for 5 × 10 5 time steps in the production stage. During this stage, snapshots are stored every 500 time steps, a total of 1000 snapshots, for further analysis and estimation of the sessile drop contact angle value. In the case of water on the LJ 9-3 surface, 4000 water molecules are considered in a simulation box with dimensions 112 Å × 112 Å × 81 Å, except for the case ε sf = 1.8 kcal/mol for which a larger simulation box with dimensions 180 Å × 180 Å × 81 Å was required to stabilize the sessile drop. For the water sessile drops, simulation times are of the order of ∼70-140 ns. The reasons for these long simulation times are explained in detail in Sec. V C. For the system of water on graphite, 2197 water molecules are placed in a simulation box with dimensions of 88.56 Å × 89.48 Å × 90.00 Å. These lateral dimensions correspond to 6048 carbon atoms in total, arranged in two graphene layers separated by 3.42 Å in the z-direction.

A. Identification of the LV interface
Instead of the traditional number-density calculation, carried out by slicing the simulation box into bins and counting the number of molecules in each bin, Willard and Chandler 43 have proposed the use of Gaussian convolution for spatial coarse-graining to calculate the (coarse-grained) density and identify the interface. The spirit of their method is similar to the smoothed particle hydrodynamics (SPH) approach applied at the continuum level, proposed by Gingold and Monaghan 54 and used by a number of other authors. 55,56 SPH renders a smooth and continuous field from a given set of discrete points; for example, from a set of atomic coordinates, it provides a smooth and continuous density field. Instead of the Gaussian kernel used in this study to obtain coarse-grained density, other kernels could also be used. 57,58 Willard and Chandler used the coarse-grained density obtained with their method to identify the water LV interface and also the liquid-liquid interfaces near melittin dimers. Here, we adopt the same approach to calculate the coarse-grained density at any point in the simulation box. In addition, the Willard and Chandler approach provides a straightforward way to calculate the outward normal to the average LV interface, a key quantity in our methodology to estimate the sessile drop contact angle. Furthermore, the Willard and Chandler approach is easy to implement and provides a continuous density field. Though this approach avoids the problem of selecting an appropriate bin size, an appropriate choice of coarse graining length is needed and it is chosen based on the far-field density profile convergence. Specifically, the coarse-grained density at position r and time t is given bȳ where the sum is taken over the number of fluid molecules N, φ(r; ξ) is the choice of spatial coarse-graining, and ξ is the coarse-graining length. A convolution with a normalized Gaussian function given by 43 is chosen for the spatial coarse-graining. The density at any position at a given time is computed using Eq. (2), and the iso-density contour with the average density ρ (r) equal half of the bulk liquid density is identified as the average LV interface. Although the summation in Eq. (2) is over all molecules, the contribution of molecules at larger r is negligible and the summation is limited to particles satisfying the criterion r ≤ 3ξ, to reduce the computational cost. 43

On the choice of ξ
Computation of the coarse-grained density using Eqs. (2) and (3) requires an appropriate choice of ξ. In this study, the mean-density profile of the fluid with respect to the instantaneous LV interface (the planar LV interface on average) is used for the selection of ξ. The simulation setup for the meandensity profile calculations is the same as that used for the LV interfacial tension calculations in Sec. IV, i.e., liquid and vapor phases coexisting with two planar LV interfaces. For a certain trajectory in the simulation and a given ξ value, an instantaneous LV interface is identified as the two-dimensional manifold with coarse-grained density [Eq. (2)] as half of the bulk liquid density. The proximity of a fluid molecule to the instantaneous LV interface, at time t, is computed using 43 where s(t) represents a point on the instantaneous LV interface identified, r i (t) is the position of the ith molecule, s * (t) is the nearest point to r i (t) on the instantaneous LV interface, andn(t) is the unit vector corresponding to the normal n(t) = −∇ρ(r, t)| r=s(t) on the instantaneous LV interface.
The mean-density profile as a function of normal distance from the instantaneous LV interface (denoted as d) is calculated as 43 where L is the length of the simulation box in the direction parallel to the LV interface. For different ξ values, ρ inst (d) is calculated and a ξ value is chosen for which the far-field value of ρ inst (d) does not change with further increase in the value of ξ (and is also found to be close to the bulk liquid density). This ρ inst (d) is obtained by averaging over 10 000 snapshots in the production stage.

B. Calculation of the contact angle
Once the average LV interface of the sessile drop is identified, two methods are used to compute the sessile drop contact angles. The first method (Method I) is the usual spherical-cap approximation and the second method (Method II) is more general in that it involves the calculation of the angle between the (outward) unit normal to the average LV interface and the unit normal to the solid surface. Clearly, Method II does not require the a priori spherical-cap approximation of the LV interface.

Method I
In the first method, the average LV interface of the sessile drop is assumed to be a spherical cap and is fitted, using least squares regression, to the equation of a sphere, i.e., where (x c , y c , z c ) and r are the center and radius of the sphere, respectively, fitted to the average LV interface and (x, y, z) represent points on this interface. A circular contour [cf. Fig. 1(a)] is chosen from the sphere centered at (x c , y c , z c ), and the contact angle (θ I ) is calculated using tan α = z c √ r 2 −z 2 c and θ I = α + 90 • . As mentioned previously, simulation data closer to the solid-surface than a cutoff of r c = 2.5 σ ff are ignored to avoid the influence of large density fluctuations.

Method II
In the second method [cf. Fig. 1(b)], the angles formed by the (outward) unit vector normal to the average LV interface (ˆ n ) with the unit vector normal to the solid surface (n solid surface , parallel to the positive z-axis) are calculated as where the unit normal is given as The outside average in Eq. (7) is taken over interface points and as a function of normal distance from the solid surface (z), and θ II = β(z = 0) is the contact angle. The presence of large density oscillations close to the wall limits the direct use of Eq. (7) at the three-phase contact line (z = 0) to estimate the contact angle. Instead, molecular simulation data close to the solid surface are ignored, and the angle in Eq. (7) as a function of z is computed, fitted, and extrapolated to the three-phase contact line to estimate the sessile drop contact angle (similar to the corresponding calculation in Method I). But in Method II, we have not made any assumptions regarding the shape of the average LV interface. It is noteworthy that the linearity of cos β(z) vs z would reflect a constant curvature configuration, hence a spherical cap. A circular contour is chosen from the sphere of radius r centered at (x c , z c ). tan α = z c / r 2 − z 2 c and θ I = α + 90 • ; (b) Method II. β = arccos ˆ n ·n solid surface is the angle between the unit vector normal to the average LV interface ( n ) and unit normal to the solid surface (n solid surface ). The filled black circles represent points on the LV interface, and the densely dotted line represents a part of the interface.

IV. TEST-AREA SIMULATIONS
For comparison, the LV, solid-liquid, and solid-vapor interfacial tension values required in Eq. (1) are obtained using the test-area perturbation method. 10,12 For a system with N particles in a volume V and at a temperature T, the interfacial tension γ can be obtained from the change in free energy as a result of change in the interfacial area A of the corresponding interface at constant volume, i.e., The change in the Helmholtz free energy can be calculated as 10 where ∆U = U 1 U 0 is the change in the internal energy of The system is equilibrated for 20 × 10 5 time steps followed by 100 × 10 5 time steps in the production stage, with a time step of 2 fs. At every 100th time step during the production stage, the energy change as a result of perturbation in the interfacial area is calculated and used in Eq. (10) to compute the interfacial tension values.

A. Test-area perturbation method
We first validate our test-area perturbation calculations by comparing our computed LV interfacial tension values with those in the literature. In the case of the LJ fluid, we find that our computed LV reduced interfacial tension, γ * (= γσ 2 ff /ε ff ) = 1.06 ± 0.01, is in good agreement with the values obtained by Gloor et al. 10 (γ * = 1.086 ± 0.023) and Peng, Birkett, and Nguyen 21 (γ * = 1.089 ± 0.018) at the same reduced temperature T * = 0.7. For SPC/E water, at T = 300 K, our computed LV interfacial tension γ = 59.91 ± 1.10 mJ/m 2 is also in good agreement with the value of 60.8 mJ/m 2 obtained by Vega and de Miguel. 46 The same method is used to compute the solid-liquid interfacial tension for the LJ and water systems and the solid-vapor interfacial tension for the water system. For water on LJ 9-3 surface, two additional test-area perturbation simulations are performed to obtain γ sv + 2γ lv and γ sl + γ lv . These two quantities in conjunction with γ lv are used to calculate the Young's contact angle (θ Y ). Based on the previous study by Peng, Birkett, and Nguyen, 21 the solid-vapor interfacial tension of the LJ system is assumed to vanish, as its contribution is negligible at T * = 0.7. We have also considered different substrate-fluid interaction strengths (ε sf ) in order to model different contact angle scenarios. The variation of θ Y as a function of ε sf is shown in Figs. 2(a) and 2(b) for both the LJ fluid and SPC/E water systems, respectively, in contact with a LJ 9-3 surface. The corresponding contact angle values can be found in Tables I and II.

B. LV interface of the sessile drop
Computation of the coarse-grained density and identification of the sessile drop LV interface requires a ξ value in Eqs. (2) and (3). An appropriate value is chosen with the help of the same setup as that used in the test area method for the calculation of the LV interfacial tension, i.e., a stabilized liquid block in contact with two vapor phases forming two planar LV interfaces. The mean-density profile with respect to the instantaneous interface (ρ inst ) is calculated, following the procedure  ) for a LJ fluid on the LJ 9-3 surface at the reduced temperature T * (=k B T /ε ff ) = 0.7. The system contains 5072 LJ particles. r * I and θ I are the radius of the fitted sphere and contact angle computed using Method I, respectively, and θ II is the contact angle computed using Method II. ε * sf (= ε sf /ε ff ) is the reduced fluid-solid potential well depth.  ) for the SPC/E water on a LJ 9-3 surface at T = 300 K. The system contains 4000 water molecules. (θ II ) and r II are the contact angle and radius of the water sessile drops estimated using Method II. ε sf is the oxygen (of water) interaction strength with the solid-surface.  Fig. 3. In the case of the LJ fluid, oscillations in ρ inst close to the interface are evident for ξ = 0.7; these are suppressed with larger values of ξ, as can be seen in Fig. 3. The far-field value of ρ inst also changes with increasing values of ξ and approaches the bulk density value ρ * = 0.835 for ξ = 1.3 at T * = 0.7. A further increase in the value of ξ to 1.5 does not result in a significant change in the far-field value of ρ inst . Thus, the value of ξ = 1.3 is chosen to compute the coarse-grained density at the same temperature T * = 0.7. In the case of SPC/E water sessile drops, with bulk liquid density value of 0.033 molecules/Å 3 at T = 300 K, ξ = 2.6 is chosen. For SPC/E water, ρ inst in Fig. 3 are very similar to those obtained by Willard and Chandler. 43 Iso-density contours with average density ρ * ≈ 0.420 at T * = 0.7 for the LJ fluid and ρ ≈ 0.017 molecules/Å 3 at T = 300 K for SPC/E water are identified as the average LV interfaces and are adopted for further analysis to estimate the sessile drop contact angle values.

C. Contact angle and relaxation dynamics
We first consider 5072 LJ fluid particles in a simulation box with dimensions 60 × 60 × 50 (in units of σ ff ) on a LJ 9-3 substrate at T * = 0.7 for the sessile drop simulations. The average LV interface of the LJ sessile drop is identified (as discussed in Sec. III A) and fitted to a sphere (Method I, as described in Sec. III B) to estimate the sessile drop contact angle values. These values are in good agreement with Young's contact angle values (Table I), as also observed in previous studies. 18,20 A snapshot of the LJ sessile drop on the LJ 9-3 surface is depicted in Fig. 4 and the spherical nature of the interface is clearly seen.
We then proceed to apply Method II (Sec. III B) to obtain the cos β vs z profile and extrapolate it to the three-phase contact line at z = 0. The sessile drop contact angle values found from Method II are also shown in Table I and also found to be in good agreement with Young's contact angles. The cos β vs z data from Method II (Fig. 5) can be seen to follow a linear profile for all contact angle scenarios at T * = 0.7. Moreover, the similar cos α vs z behavior from the spherical fit (according to Method I) to the LV interface is also provided in Fig. 5 for comparison. Summarizing, Method II is able to capture the spherical-cap nature of the LV interface, connected to a linear behavior of cos β vs z data, without the need for a priori assumptions on the shape of the average LV interface.
Having confirmed the validity of our proposed method for sessile drop contact angle values, in the case of the LJ fluid on the LJ 9-3 surface, we now consider SPC/E water on the LJ 9-3 surface at 300 K for the sessile drop simulation (total 4000 water molecules setup in a lattice configuration, initially). The water sessile drops took 70-100 ns to relax to a spherical cap, depending on the fluid-surface interaction strength. The LJ sessile drops on the LJ 9-3 surface were also initiated in a lattice configuration, but only took 10-15 ns to reach the spherical shape. It should be noted that the time scale observed for the water sessile drops to become fully spherical, 70-100 ns, is quite different to the time scale of 20-30 ns reported in the literature. [27][28][29][30][31][32][33] This could be most likely due to the chosen initial configuration, a lattice configuration. Another common initial configuration is a pre-equilibrated spherical drop which is then brought in contact with the solid surface; such a procedure is likely to result in a faster relaxation of the sessile drop to the spherical shape. For water drops to equilibrate to spherical shape on a graphite surface, Shih et al. 41 reported a time scale similar to that in our study, ≈100 ns, but the initial configuration of the water molecules was not mentioned in their study. As already emphasised, a linear trend of the cos β vs z profile (cf. Fig. 6) of equilibrated water sessile drops reflects the spherical nature of the LV interface. Furthermore, the estimated water sessile drop contact angles are in good agreement with Young (cf. Table II). These results suggest that sessile drops with a number of molecules of the order of thousand can be used to extract good estimates of the contact angle. The use of a LJ fluid on a LJ 9-3 surface and the testarea method has enabled us to verify the validity of Young's equation. A number of studies [20][21][22]25,26 have followed this approach, but only a limited number of them have done the same analysis for water sessile drops. Furthermore, the comparison between the sessile drop contact angle values and Young's values has helped us shed light on the ambiguity that exists in terms of system size and simulation time for water. It should also be noted that in Tables I and II we report the error estimates for Young's contact angle, but not for the sessile drop contact angle values. This is because the calculation of the latter is influenced by a large number of parameters such as bin size or coarse-graining length (ξ) involved in the local fluid density calculation, a chosen density value of the iso-density contour to identify the LV interface, number of configurations chosen for the sessile drop analysis, and the cut-off distance chosen to ignore data close to the substrate. The error in the radius of the sessile drop in Tables I and II only represents the error from fitting the LV interface to a sphere and does not include the effect of the parameters mentioned above.
Finally, in the case of an SPC/E water sessile drop at 300 K on a graphite surface, a time of ≈140 ns was required for the drop to attain spherical shape, with 2197 water molecules initially setup in a lattice. Four snapshots at different times during the simulation are shown in Fig. 7. In addition to the final linear trend of the cos β vs z data at 136 ns, which confirms the spherical-cap nature of the interface for large times, the cos β vs z data at intermediate simulation times are also shown in Fig. 8. The evolution of the sessile drop toward the spherical shape can be seen in Fig. 7 as also reflected in Fig. 8 where the transition of cos β vs z profile from non-linear to linear is evident. A number of intermediate contact angle values during the drop evolution to the spherical shape can also be estimated using our method. For these intermediate sessile drop contact angle values, a partial set of cos β vsz data (for which a linear fit can be assumed) is fitted to a straight line and the fit is extrapolated to the three-phase contact line (z = 0). The intermediate sessile drop contact angle values, thus obtained, are in good agreement with the values estimated from the spherical sessile drop at 136 ns (cf . Table III).
In a previous study, to understand the effect of surface polarity on water, Giovambattista, Debenedetti, and Rossky 35 fitted the r(z) vs z profile of SPC/E water drops to a generic FIG. 7. Snapshots of water drops on graphite (composed of two graphene layers separated by 3.42 Å in zdirection) at different times with simulation temperature T = 300 K. Initially, water molecules are set up in a lattice configuration on graphite. Carbon atoms of the graphite surface are kept fixed during the simulation. (a)-(d) correspond to the simulation times of 37, 72, 100, and 136 ns, respectively. quadratic fit, considering only data below r(z) > 30 Å. The non-sphericity of the water drop may have been the reason for the contact angle estimation with only limited data, although unfortunately the finer details were not mentioned by the authors. In our study, a partial set of cos β vs z data close to the substrate (that follows linear trend) is a reasonable choice and results in good estimates of the sessile drop contact angle, as shown. The extrapolated value of cos β vs z data at the three-phase contact line (z = 0) at intermediate times (on-route to the spherical shape of the interface) intersect close  to the same point (cf. Fig. 8), revealing a constant-contactangle relaxation dynamics. This finding is also consistent with the results by Lukyanov and Likhtman 59 that interfacial tensions have relaxation times faster than those for the drop shape. (The idea of interfacial-tension relaxation is at the heart of a macroscopic hydrodynamic model introduced by Shikhmurzaev 60 and was scrutinized asymptotically by Sibley, Savva, and Kalliadasis. 61 ) Although not reported here, we observed similar constant-contact-angle relaxation dynamics for the water sessile drops on the LJ 9-3 surface. Upon this observation and with our proposed methodology, one does not need to wait for a water drop to reach the spherical shape to achieve good estimates of the contact angle; short simulation times such as 10-20 ns are sufficient (in Table III). Furthermore, the maximum number of water molecules considered is only 4000 molecules and was still able to obtain contact angles in good agreement to Young's equation (Table II). This is a drastic improvement, for instance, to the results of Santiso, Herdes, and Müller 34 who concluded the requirement of half a million water molecules to obtain good agreement with Young. Moreover, the spread in the angle distribution obtained from our methodology is small (cf. Figs. 5, 6, and 8), as expected for a homogeneous solid surface, compared to the spread obtained by Khalkhali et al. 39 Our method and those of Santiso, Herdes, and Müller 34 and Khalkhali et al. 39 are somewhat similar in spirit, in that the spherical-cap approximation is not invoked, instead the angle between the (outward) unit normal to the average interface and the normal to the solid surface is used; but our method is shown to be robust as well as easy to implement.

VI. CONCLUSIONS
A new methodology is proposed for the direct measure of sessile drop contact angles following the Willard and Chandler 43 method to calculate the coarse-grained density and outward normal to the average LV interface. The methodology is exemplified by calculating the contact angle of a LJ fluid drop on a LJ 9-3 surface and of an SPC/E water drop on LJ 9-3 and graphite surfaces. In all cases, we find good agreement with Young's equations. The interfacial tensions required in Young's equation are obtained by employing testarea simulations for the corresponding planar interfaces. We have been able to resolve a number of previous issues that still eluded us such as the ambiguity in terms of system size and simulation time. In particular, our study demonstrates that it is neither the system size nor the simulation time that limits the accurate estimation of the contact angle, but rather the particular method adopted for its calculation. For instance, the use of the spherical-cap approximation for a water drop requires long simulation times of the order of 100 ns, but this is not necessarily the case for the equilibration of the contact angle. Although the entire water sessile drop does not adopt the spherical shape even for long simulation times such as 100 ns, the contact angle is well established much faster, i.e., within 10-20 ns of simulation time, even when starting from a configuration of water molecules in a lattice. As a matter of fact, the sessile drop evolves toward the spherical shape, following constant-contact-angle relaxation dynamics.
It would be of interest to study several related problems, for example, analyzing the detailed characteristics of the dynamics of drop evolution including the drop shape, density fluctuations in the vicinity near the solid-fluid interface, and the layering effect on the bulk fluid, but also to extend our methodology to other directions, such as variation of contact angle with drop size, and other physical settings such as dynamic wetting on planar 62,63 and topographically and/or chemically heterogeneous substrates. 62,[64][65][66][67] For dynamic wetting, in particular, the spherical-cap approximation is crude and would be interesting to apply the proposed methodology for the variation of contact angle with time. For wetting of heterogeneous substrates, the methodology proposed here would allow us to obtain a distribution of the contact angle as a function of the location on the substrate rather than a single contact angle. Of interest would also be extensions of this study to LV interfaces in confinement, such as that between two parallel walls (e.g., Ref. 68). We shall examine these and related issues in future studies.