On the Calculation of Solid-Fluid Contact Angles from Molecular Dynamics

A methodology for the determination of the solid-fluid contact angle, to be employed within molecular dynamics (MD) simulations, is developed and systematically applied. The calculation of the contact angle of a fluid drop on a given surface, averaged over an equilibrated MD trajectory, is divided in three main steps: (i) the determination of the fluid molecules that constitute the interface, (ii) the treatment of the interfacial molecules as a point cloud data set to define a geometric surface, using surface meshing techniques to compute the surface normals from the mesh, (iii) the collection and averaging of the interface normals collected from the post-processing of the MD trajectory. The average vector thus found is used to calculate the Cassie contact angle (i.e., the arccosine of the averaged normal z-component). As an example we explore the effect of the size of a drop of water on the observed solid-fluid contact angle. A single coarse- grained bead representing two water molecules and parameterized using the SAFT-γ Mie equation of state (EoS) is employed, meanwhile the solid surfaces are mimicked using integrated potentials. The contact angle is seen to be a strong function of the system size for small nano-droplets. The thermodynamic limit, corresponding to the infinite size (macroscopic) drop is only truly recovered when using an excess of half a million water coarse-grained beads and/or a drop radius of over 26 nm.


Introduction
Wetting is the ability of a liquid to maintain contact with a solid surface, resulting from intermolecular interactions when the two are brought together. The presence of a liquid drop on a rigid surface is a reflection of the force balance between adhesive and cohesive forces and is commonly used to determine the wettability (the degree of wetting) of the solid-fluid system in terms of the solid-fluid contact angle, θ (see Figure 1). In this context, hydrophobicity is commonly referred to as the ability of a solid surface to repel water: if the water contact angle is smaller than 90°, the solid surface is considered hydrophilic and if the water contact angle is larger than 90°, the solid surface is considered hydrophobic. Despite the fact of being such a well-defined problem the amount of conflicting (both theoretical and experimental) reported values for a given system is intriguing (see Figure 2, data from reference [1]). For instance, in the case of the graphite-water system, contact angles have been addressed extensively by experimental and theoretical approaches; however, a single general value has not been accepted [2][3][4][5][6][7][8][9].
A variety of causes for the discrepancies can be enumerated: heterogeneity and/or impurities at the surfaces or in the fluids, different methodological unstandardized calibration of equipment in experiments, possible system size effects and the distinct interaction potentials in simulations.
At the molecular scale, the main hurdle is that estimating contact angles for nanodroplets on surfaces is complicated by the fact that there are significant fluctuations in the shape of the droplet, and its geometry at a given step is often not axially symmetric. Furthermore, for very small nanoclusters, the fluid interfacial tension is a function of the curvature, and the planar limit is, in some cases, not recovered even after drop radii of 14 times the molecular diameter [10,11]. The change in the line tension with curvature (discussed in the latter part of this manuscript) is also an important factor affecting the result. It is seen that the contact angle will be, in the nanoscopic limit, a strong function of the system size [2]. In the analysis of simulations, contact angles [12,13] are commonly determined by using two-dimensional slices of the droplet and fitting its density profile to an empirical function, usually a circular section [14]. Such an approach, although appealing from the simplicity of the method, provides inconsistent results, particularly for small droplets [11]. Understandingly, different methods of increasing complexity have been devised for this purpose [15][16][17]. Figure 3 illustrates the difficulties associated with defining the contact angle using two-dimensional slices of molecular snapshots, especially for small droplets. Such droplets are often asymmetrical, and their shape changes substantially with time, due to capillary fluctuations. Using one or a few two-dimensional projections of a few molecular snapshots can thus potentially lead to large errors on the measured contact angle: for example, the values shown on Figure 3 vary over range of about 10°.
In addition to the inherent fluctuations of the drop's shape, there is a degree of arbitrariness in defining a single line separating the two phases in the two-dimensional projection. The interface layer is actually diffuse, as seen in Figure 3, spanning a width of several molecular diameters. This introduces another potential source of error, as the contact angle measured can change over a range of 5° or more depending on the choice of the contact line.
In this paper, we propose a method that does not make any a priori assumption about the shape of the drop, and uses the complete three-dimensional structure of the droplet near the surface to estimate the contact angle. For such analysis, we propose a geometrical estimation of the contact angle based on cloud point data breakdown of an equilibrated molecular dynamics (MD) trajectory of a drop on a given surface, this eliminates the effect of both spatial and temporal fluctuations on the estimated value of contact angle, yielding a value that better represents the average shape of the drop, as explained in detail in the next section.

Methodology
In order to estimate the average contact angle during an MD simulation, we first define a contact layer at each time step. This layer is defined as the set of molecules within the liquid-vapor interface that are close to (i.e., within a maximum distance z max from) the solid surface. We then estimate the normal to the interface for each molecule in the contact layer by finding the plane that best fits the local shape of the interface at that point. The average of all the local normal vectors for all the time steps, , can be used to calculate the average contact angle as , where is the normal to the solid surface.
In practice, at each time step in the simulation, we carry out a two-step calculation: In the first step we use a discretized density profile to identify the molecules belonging to the liquid-vapor interface. In the second step, we estimate the local surface normal vectors for the molecules in the contact layer. The procedure is explained in detail below.

Identification of Interfacial Molecules
We identify the molecules belonging to the liquid-vapor interface using the following procedure: (1) We divide the simulation box into subcells, and calculate the local number density in each subcell.
(2) We mark subcells as liquid if their number density is greater than a cutoff value , otherwise we mark them as vapor cells. The cutoff density is chosen close to the average between the liquid and vapor number densities at the conditions of interest. (3) Every subcell that is adjacent to at least one liquid cell and one vapor cell is marked as an interface cell. (4) Finally, all molecules contained within the above-determined interface cells are marked as interfacial molecules.
The results obtained using the interface-sensing procedure are only weakly dependent on the choice of cutoff density, as long as this density is between the liquid and vapor densities. The number of subcells, on the other hand, should be chosen carefully: having more subcells identifies the interface with a higher resolution. However, if the number of subcells is too large, the density fluctuations within the droplet may cause the algorithm to incorrectly label too many cells as part of the interface. This can be easily detected by visualizing the interface molecules at one or a few simulation steps. Alternatively, one can use the average coordination number or the Minkowski-Bouligand dimension (obtained from a box-counting algorithm) [18,19] as a quantitative measure to select the appropriate subcell size. As the cells become smaller and the interface-sensing algorithm starts fitting the molecular-level fluctuations in the density rather than the actual interface, both the fractal dimension and the average coordination number will start deviating from the value expected for the larger subcells, indicating that the subcell size has become too small.

Estimation of Local Contact Angles
Having identified the molecules belonging to the interface, we choose a subset of the interface molecules as the interface contact layer. This layer contains the interface molecules that are within a given distance z max from the solid surface. We then estimate the local contact angle at the position of each molecule i using the following procedure: where denotes the outer (Kronecker) product. (e). Find the eigenvector of Ω corresponding to its smallest eigenvalue. This is the normal to the plane that best fits the set of molecules in step (c). The sign of this normal is chosen to point away from the center of the droplet [20].
We average all the local normal vectors obtained at all the time steps to obtain the average normal. The component of this average normal perpendicular to the solid surface is the cosine of the average contact angle. The eigenvalues of the covariance matrix measure the variability of the positions along the three orthogonal directions defined by the eigenvectors. These can be used as a measure of how 1-, 2-or 3-dimensional the interface is. If one of the eigenvectors is close to zero, the interface is close to planar in that region.
The procedure described above to estimate the local contact angle requires arbitrarily setting a distance from the solid surface, z max , to define the contact layer, and a cutoff radius to estimate the tangent plane. In theory the contact angle would correspond to the limit when z max → 0. However, if z max is chosen to be too small, the local density fluctuations within the contact layer cause the distribution of local contact angle values to be too noisy. Thus, one should choose the smallest value that gives a reasonable estimation error for the average contact angle. For the examples shown in this work, we have found z max = 0.5 nm to give reasonable results.
For the cutoff radius, r c , a compromise is also needed: if the value is too small, local density fluctuations cause the error in the estimated normal to be too large, whereas if the value is too large the orientation of the tangent plane will be affected by molecules far from the position of interest. The choice of r c should be made by considerations similar to those used in choosing z max . For the examples in this work, we have found r c = 1 nm to be a reasonable value. A more detailed discussion of how to choose this cutoff can be found in reference [21]. Finally, we estimate the statistical error associated with the average contact angle computed by our method, we use the bootstrapping method [22].

Molecular Dynamics Details
Atomistic MD simulations are generally limited to nanometer-sized water droplets [2] Consequently, the apparent contact angle is usually drop-size dependent. To explore bigger systems, and aiming to find the optimal size for the MD contact angle calculation, we adopted a coarse-grained (CG) approach to describe the solid-water system in order to test the system size dependence beyond the atomistic limits, in a reasonable time [2].
We used the GROMACS simulation open source [23] suite to calculate the MD, which is well suited to implement Mie potentials [24]. Here, a single CG isotropic bead represents two water molecules [25]. Although several options are available for choosing the number of water molecules in a CG bead [26], this choice and the current parameterization produces sensible results, including a melting point, the surface tension, liquid densities and vapour pressures close to the experimental. The parameterization was carried out using SAFT-γ Mie approach, [27][28][29] where the water parameters were obtained by fitting to macroscopic properties, namely, the planar limit interfacial tension and liquid state density of water in a range from 0 °C to 40 °C. The SAFT EoS is a perturbation approach based on a welldefined Hamiltonian; here the CG beads are represented in the theory by a Mie potential, u: where r is the intermolecular distance, and ε, and σ, are the adjustable parameters relating to the energy and distance scales, λ a is the dispersion exponent and λ r is the short-range repulsion. The potential is expressed in terms of two constants A and C for ease of tabulating in MD codes. Solid walls are modeled implicitly and described by an integrated potential of the same form in the z-dimension. Eight implicit solid surfaces of increasing ε Wall#-W /k B increments of 10 K, from 60 K to 130 K (labeled Wall01 to Wall08 respectively) are employed [30]. Table 1 summarizes the selected coarse-grained parameters. The systems are run under a canonical (NVT) ensemble, where the total volume, concentration and temperature are kept constant. Periodic boundary conditions were applied in xy dimensions, meanwhile an attractive wall was placed at z = 0 and a repulsive one (C = 5.73512 × 10 −4 kJ mol −1 nm 4 , A = 1.93147 × 10 −6 kJ mol −1 nm 10 at the maximum height of the box. The number density of the atoms for each wall was set in 5 nm −2 (c.f. in a Si crystal the number of atoms per nm 2 on the (100), (111) and (110) planes are 6.78, 7.83 and 9.59 nm −2 , respectively [31]). The simulations are thermostated to 298.15 K every 1 ps by a Nose-Hoover algorithm, all non-bonded interactions were truncated at 2.0 nm. The trajectories were recorded every 1000 time-steps (Δt = 0.01 ps) for at least 2000 ps after equilibrium.

Results and Discussions
We simulate 16,000, 32,000, 64,000, 128,000, 256,000, 512,000 and 1,024,000 water molecules (8,000, 16,000, 32,000, 64,000, 128,000 and 512,000 beads) on the eight solid surfaces. The simulation boxes dimensions can be seen in Table 2. Although arbitrary, the reason guiding the choice of the simulation box size is to prevent the interaction between the sample and its periodic images. The meshing was done by dividing the simulation domain with a 1.2 × 1.2 × 1.2 nm 3 subcell. The water droplet density contour is obtained from the cloud point data set analysis explained in the methodology section.
To test the size dependence of the water contact angle on a given surface we chose the intermediate Wall05 (see Table 1). The contact angles obtained using a moderately hydrophilic substrate (Wall05) as a function of the drop size can be seen in the Figure 5. The water contact angle over Wall05 exhibits a marked system-size dependence even up to 256,000 water molecules (128,000 beads) Larger drops exhibit a less pronounced but nevertheless noticeable size dependence. It is interesting to note that the effect of this scale up is an increase in the hydrophobicity of this surface in correspondence with previous simulations [32] and experimental studies at the macroscale. For this substrate, a system with 1,024,000 water molecules shows an apparent limiting contact angle of 74.30° [8]. The Wall05 parameters are chosen to loosely relate to graphene, a hydrophilic substrate. Following from the above, we use the system of 256,000 water molecules to study its contact angle as a function of the fluid-substrate interactions; the results can be seen in Figure 6. As expected, increasing the interaction energy between water molecules and the attractive wall diminishes the solidfluid contact angle. For this particular coarse-grained surface model the functionality is linear, and a value of ε Wall-W /k B ~ 85 K can be taken as the boundary between hydrophilic and hydrophobic surface behavior.
An ancillary quantity of interest in the context of the drops on surfaces is the line tension. The line tension is the relation between the energy associated with the three phase contact line and the length of this line. This quantity is inherently dependent on the size of the drop [33] and is discussed as one of the reasons why the calculated tensions appear to be size-dependent. For large drops, one can obtain an estimate of the line tension using the approximation described by Weijs et al. [34]. For a typical case as shown above: the case of a drop of roughly 26 nm in diameter and considering the planar limit of the fluid-vapor tension (72 mN·m −1 ) results in a contact tension line strength of O(10 −8 ) J·m −1 . While this number seems to be different from those estimated from micrometer drop experiments [35] or from molecular simulation of atomistic water models [36,37]; it is appropriate to point out that values as low as ×10 −11 J·m −1 and as high as ×10 −5 J·m −1 have been reported [33], which place our results in the correct context. Obviously, there is scope for much more detailed research into this topic.

Conclusions
We have proposed and validated a methodology for the unambiguous calculation of the solid-fluid contact angle from molecular dynamics simulations. We have tested model coarse-grained water-solid systems far beyond the limits commonly taken in atomistic simulation and showed, that for this particular model, more than 500,000 effective beads and/or drop diameters in excess of 50 nm would be required in order to obtain a result which is invariant of system size. So far the methodology has been applied over homogenous surfaces. However, is a well-known fact that surface roughness and energetic heterogeneities will have a profound effect on the contact angle calculations. The methodology presented is well suited to capture those effects.