1. Introduction

“Wetting” is the ability of a liquid to maintain contact with a solid surface through intermolecular interactions. In other words, wetting is the balance between adhesive and cohesive forces; cohesive forces within a liquid cause it to avoid contact with the surface by forming a spherical shape when all three phases are in equilibrium, whereas adhesive forces between the liquid and the solid surface cause the liquid to spread across the surface. If a surface of a material is attracted to a certain liquid, it is considered to be “-philic” for that substance; for example, a surface that is attracted to water is called “hydrophilic”, and attracted to oil is called “oleophilic”. However, if a surface tends to repel a certain liquid, it is considered to be “-phobic” for that substance; for example, “hydrophobic” and “oleophobic” mean that the surface repels water or oil, respectively.

Thus, the goal of this report is to survey the use of molecular dynamics (MD) simulations for investigating the wetting behavior of liquids on solid surfaces at the nanoscale by calculating both the surface tensions of solids and liquids and contact angles of liquid nanodroplets on surfaces with molecular scale roughness. This work is first put into the context of what is done experimentally via wetting models for solid surfaces interacting with liquid droplets.

2. Brief Overview of Wetting Models for Solid Surfaces

The wettability of a surface is determined experimentally by placing a droplet of liquid on the surface and measuring the apparent contact angle (θ), which is defined as the angle formed between the solid/liquid interface and the tangent line of the vapor/liquid interface [1] and is depicted in Figure 1. Assuming that the liquid is water, if the apparent θ is smaller than 90˚, the surface is considered to have high wettability, or hydrophilic; larger than 90˚ is considered low wettability, or hydrophobic; and over 150˚, the surface is classified as superhydrophobic [2, 3]. These terms are defined similarly for other types of liquids.

Figure 1. Wetting characteristics of a solid surface, from low wettability (left) to high wettability (right).

Figure options

When the solid surface is perfectly flat, rigid, smooth and chemically homogenous (in other words, an ideal surface), the contact angle obtained from the droplet is called Young’s contact angle, which is related to the surface tension (g) of each substance through the following formula [4],

where g_SV is the surface tension between the solid and gas (vapor), g_SL is between solid and liquid, and g_LV is between liquid and gas. In general, θ increases when γ_LV decreases; in other words, the lower γ_LV is, the easier it is for that liquid to wet the surface [5]. Thus, lowering γ_SV will be useful in order to obtain better hydrophobicity and even superhydrophobicity; however there is a limitation to the lowest g that a material can possess. For organic liquids with a much lower g, only decreasing g of the surface has been observed to not be enough to achieve superoleophobicity [2, 6]. Note that g for a substance is comprised of all of the different intermolecular interactions (cohesive forces) that exist within a substance, according to the Fowkes equation [7, 8],

where γ^D is the contribution due to dispersion (van der Waals) forces, γ^P is permanent dipoles, γ^ind is induced dipoles (polarizability), γ^H is hydrogen bonding, and γ^m is metallic interactions.

Real surfaces do not exhibit perfect smoothness, rigidity, or chemical homogeneity, and thus other models have been developed to take these effects into account. One consequence is contact angle hysteresis (H) [1, 3] that occurs because there are a range of thermodynamically stable θ values (called metastable states) that may exist. H is thus defined as the difference between the maximum stable angle, called the advancing contact angle (q_a), and the minimum stable angle, or the receding contact angle (q_r), which are commonly measured on a sloped plane. In addition, surface roughness plays a significant role [1, 9, 10]. For example, the lotus leaf exhibits superhydrophobicity, and scanning electron microscopy (SEM) images revealed that a lotus leaf contains not only primary protuberances at the microscale but also each protuberance is comprised of secondary rough structures on its surface [11]. In terms of surface roughness, there are two types of wetting regimes. One regime is the heterogeneous wetting regime, where a gas (often air) fills the spaces created by the rough texture rather than the liquid, and thus the liquid-solid contact line also contains liquid-gas and solid-gas contacts. This regime is described by the Cassie-Baxter model, given by the following, [9]

wherer is the ratio of the true area of the solid surface that is wet to the apparent contact area if it were perfectly smooth, and f is the fraction of the solid surface area that is wet by the liquid. The other regime is the homogeneous wetting regime, where the liquid completely fills the spaces created by the rough texture of the surface. This regime is described by the Wenzel model, which is obtained from Eq. 3 when f = 1 [12]. When θ is fixed, θ increases monotonously with f. In other words, the higher that f is, the lower that r will be. Thus, to find an optimized r is a major goal to design superhydrophobic and superoleophobic surfaces. Also, for the Wenzel model, the contact angle hysteresis will be very high, meaning that q_r will be very low compared to q_a [10], resulti ng in difficulties when removing a liquid from the surfaces when compared with the Cassie Baxter model, which performs much better in this aspect.

Experimentally, g of the surface can be tuned through surface modifications. For example, cotton was observed to achieve superhydrophobic and self-cleaning properties with treatment with a fluorosilane [13]. Treatment with fluorinated groups is effective for lowering g because fluorine has a small atomic radius and the largest electronegativity [14, 15, 16, 17]. Another widely used approach is to enhance the surface roughness with a fractal structure to improve the hydrophobicity of surfaces [2, 18, 19, 20]. Recently, Lee and coworkers designed and developed new anti-icing textile fabrics through the combination of lowering g and increasi ng surface roughness [6]; the results suggested that the apparent contact angle of water was 152-158˚, which also agreed with the prediction from the Cassie Baxter model, indicating that those treatments led to the fabric exhibiting superhydrophobicity.

3. Applications of Molecular Dynamics Simulations for Characterizing Surface Wetting

Molecular dynamics (MD) simulations are comprised of the physical motions of atoms and molecules that are calculated by solving Newton’s equations of motion [21]. The forces are derived from the potential energy of the atoms within molecules that are comprised of both intramolecular and intermolecular forces. The intramolecular interactions are calculated based on the geometry of the molecule, such as bond energies, three body “bend angles”, and four body torsional angles, while the intramolecular interactions are a combination of van der Waals and electrostatic forces. By solving these equations at each time step, a time trajectory for a set of molecules can be calculated. MD simulations have provided a powerful tool for studying a variety of physical and mechanical properties, as this methodology provides a way to quantify the influence of atomic-scale interactions on larger scale properties. For example, MD simulations can predict surface characteristics such as wetting due to chemical modifications before applying the real treatments in the laboratory [22, 23, 24, 25].

3.1. Predicting the surface tensions of liquids and solid surfaces

As indicated in Figure 1 and Eq. 1, one critical interfacial property to characterize surface wetting is g. In an MD simulation, g can be calculated in several ways. One way is to calculate the excess energy between a molecular model of a surface (E_surface) and its corresponding bulk model (E_bulk), given by the following equation [26, 27, 28],

where A is the interfacial surface area of the surface. (Note that most MD simulation software outputs the energies on a per mole basis, and thus dimensional analysis may be necessary). Because the potential energy can be broken into the individual contributions due to factors such as hydrogen bonding and van der Waals interactions, MD simulations can provide additional insights into understanding g ( and thus the wettability) of a substance through the Fowkes relationship (Eq. 2), since those terms cannot be easily measured experimentally. Another way to calculate g from an MD simulation is by calculating the pressure profile using the following equation [29, 30],

where P_xx, P_yy and P_zz are the average pressures in the x, y and z dimensions, respectively. Finally, g_SL is also directly proportional to the work of adhesion (W^adh) at the solid-liquid interface through the Dupre and Dupre-Young equations defined, respectively, as [7],

In other words, W^adh describes the work required to separate the solid and liquid by creating unit areas of each (defined by g_SV and g_LV) but at the expense of the solid-liquid interface (g_SL). Thus, a large positive W^adh value signifies strong adhesion between the solid and liquid, whereas a smaller positive value suggests weak adhesion; a negative value indicates repulsion, suggesting that the liquid will spontaneously separate from the solid. From an MD simulation, W^adh can be calculated from the following equation [31],

The E_inter term, which is the interaction energy between the solid and liquid, is calculated by first building a molecular model where the surface is completely covered with a liquid, and then calculating the total potential energy of the system (E_total), as well as the total potential energy of the isolated solid (E_solid) and isolated liquid (E_liquid) components at a particular MD time step. This approach has been used to determine the relative differences for particular liquids interacting with surfaces based on different surface conditions such as patterning and roughness [32, 33, 34].

3.2. Predicting contact angles of liquid nanodroplets

MD simulations have also been used to study wetting through the use of a liquid nanodroplet in contact with a surface. This technique is similar to the sessile drop technique used in physical experiments but on a different size scale. This calculation is generally done by first building a molecular model of a droplet of the liquid, and then equilibrating that droplet while it is in contact with a molecular model of a surface of interest with a particular degree of molecular surface roughness. Once a minimum energy state is achieved, the shape of the liquid is examined; if the liquid forms some sort of a spherical (or lenticular) shape, as with physical experiments, a static θ (depicted in Figure 1) may be calculated, and this value can be correlated with g_SL via Eq. 1. As will be discussed below, there are some challenges with calculating θ from nanodroplets; however, these types of MD simulations have still provided some significant insights at the nanoscale on wetting characteristics that have complemented experiments.

One significant challenge with nanodroplets from MD simulations is that the conventional θ becomes ill-defined at the molecular level, even for smooth surfaces, as illustrated in Figure 2. However, different methods to approximately fit the outline of a droplet onto a spherical shape have been utilized [26, 35, 36, 37, 38, 39]. One approach is to neglect the liquid molecules that are closest to the surface, because there is frequently restructuring that occurs within the liquid molecules closest to the surface due to the strength of the interaction between the surface and the liquid; typically, this distance is 8 - 10 Å [40, 41, 42, 43]. This approach is done by physically drawing a line that is equal to the slope of the droplet just above the neglected thickness, then using this line to calculate the contact angle between the droplet and the surface. The advantage of this technique lies in its simplicity, as it only requires a liquid density profile and a way to draw a line. The disadvantage is the lack of an empiri cal method to verify the accuracy of the line placement, thus allowing for some error to enter into the calculation of θ. Considering the anisotropy of the wetting behavior, θ values are often calculated from different side views and averaged [25, 44].

Figure 2. Depiction of the variation in the contact angle for a nanodroplet.

Figure options

Another method that is widely used is to detect the shape of the droplet from the density profile [29, 37, 44, 45, 46]. For example, Giovambattista and coworkers [39] divided the droplet into slabs of width of 0.5 nm, separated vertically by 0.25 nm, and parallel to the surface; for each slab, they calculated the horizontal density profile. For each slab centered at the z axis, they associated a drop radius, r_drop (z), given by the distance at which the horizontal density profile falls to 0.2 g/cm³. The contact angle is obtained from the droplet profiles by fitting the curves with a quadratic function. By this method, the authors achieved smoother curves of the outline of the droplet and then measured the contact angles from it. Another method was based on a mathematical relation between the average height of the center of mass of the water droplet (<z_c.m.>) and contact angle, given by the following equation [37],

where R₀ is the radius of an isolated spherical droplet of N water molecules.

Another challenge with correlating MD simulations to experimental results is that there is no consistent methodology for the preparation of the droplet or the surface, and using different initial configurations can cause a wide disparity in the resulting values. An example of this disparity is the wide range in published values for θ of water on graphene from both experiments and simulations, which range from 29° to 115°, with the most common range reported being in the range of 84°-86° [47]. Studies have indicated that basic system parameters such as the size of the nanodroplet can have an influence on the calculated θ values. In a study of θ of water on a gold substrate, Wu and coworkers [48]determined that changing the droplet diameter from 4.5 nm to 7.5 nm resulted in a change in θ from 54° to 117°. This dependence on the droplet size makes it difficult to compare reported θ values without a consistent droplet size.

Even the formation and placement of the droplet on the surface is inconsistent among different studies. A common initial liquid configuration is to create a cubic simulation box of the liquid and , before placing the liquid in the interaction regime of the surface, allow the liquid to form a spherical shape through an MD simulation equilibration of just the liquid [39, 40, 41, 49, 50]. An alternative methodology is to equilibrate a cubic simulation cell of water at the desired temperature and pressure, and then placing this cube of water a small distance from the surface that has significantly larger dimensions [42, 51, 52, 53]; MD simulations are then applied to the combined system under the desired thermodynamic conditions for a period of time that is sufficient for enabling the cube to become a droplet interacting with the surface. Both methodologies seem to have their own advantages and disadvantages. The advantage of using the cubic liquid cell on the surface is the decrease in computational time as compared to the spherical nanodroplet on the surface. Alternatively, the advantage of using the spherical nanodroplet is that it is more representative of real world experimental conditions.

Despite these difficulties, MD simulations have provided insights at the nanoscale of what affects wetting behavior at the atomic level, like droplet size [25, 45], surface polarity [39], and surface structure [25, 26, 34]. The ability to easily make systematic small or large changes to the systems in MD simulations and to quantitatively calculate the influence of these gives an understanding of the mechanisms that govern surface wetting [25, 26, 35, 39, 54]. For example, although it is common to change the composition of the surface while keeping the nonbonded interaction parameters between the liquid and surfaced fixed, a number of studies have changed the strength of interaction between a liquid and a droplet through changing the strength of the nonbonded interaction parameters [32, 55, 56, 57]. MD simulations have also been used to examine surface properties such as surface roughness and their influence on the spreading of the nanodroplet [41, 58]. An extreme example of this modification is the formation of surfaces that are patterned with nanoscale grooves [59] or pillars [52]. The pattern can be easily modified to investigate the influence of the height and spacing of the grooves/pillars. The simulations allow for precise control of the placement of the droplet on the patterned surface, and allow for precise measurement of the extent to which the droplet penetrates the space between the groves/pillars, and thus a systematic investigation of the Wenzel versus the Cassie-Baxter models can be performed.

MD simulations have also been used to calculate the moving contact angle [40]. By applying a directional force to the liquid parallel to the surface, the droplet can be distorted similar to a droplet on an angled plane. The droplet no longer has a contact angle equal in all directions, but now has an advancing and receding contact angle. These two angles, specifically, the difference between these two angles, can be used to calculate the capillary number of the liquid with the surface; this number can be used to relate the strength of interaction between the liquid and the surface. By then removing the external force, the contact angle hysteresis may be investigated.

Finally, one adaptation of contact angle simulations with MD is its use to refine molecular mechanics force fields. For example, Jaffe and coworkers [32] used the contact angle between water and graphite to refine their existing force field parameters for carbon-water interactions. The authors began with established water (SPC/E) and carbon (continuum model) nonbonded parameters, and then they modified the cross-term nonbonded parameters between carbon and the oxygen in water to match the contact angles reported from experimental results. This example highlights the usefulness of the wettability calculations for also further improving computational results.

4. Future Outlook

As technology advances, computational resources become significantly cheaper and more advanced, thus enabling the simulation of systems that are increasingly complex with less computational time. Larger droplets and surfaces could be used, negating some of the variations in droplet surface density in systems with a small number of liquid molecules. Additionally, increasingly complex computational models that account for forces such as polarization and droplet dynamics can be applied to the system and improve the accuracy and information garnered from the simulations.

By applying a force in a certain direction in simulations to mimic roll-off droplet experiments, it would be insightful to not only characterize advancing and receding contact angles, but to also quantify the amount of residue left by liquids in the process, or the robustness of the wettability [2]. In addition, as described above, the formation of a droplet on a surface is dependent on the interplay between the adhesion of the liquid with the surface and the cohesion of the liquid with itself. Typically most studies focus on the adhesion, neglecting further insight into the cohesion of the liquid and its effects. This area could be an interesting focus of future research, especially with the use of polarizable force fields. Finally, the environment of wetting in MD simulations is usually in vacuum, while in experiments they are conducted in air, and thus investigating how the wetting may be impacted by the incorporation of explicit gas molecules would be another useful future direction, especially in the context of the wetting models described by Eqs. 1-3.

Acknowledgments

This material was sponsored by the Defense Threat Reduction Agency - Joint Science and Technology Office for Chemical and Biological Defense (grant number HDTRA1-08-1-0049, HDTRA1-12-C-0048). The U.S. government is authorized to reproduce and distribute reprints for governmental purposes notwithst and ing any copyright notation thereon. DISTRIBUTION STATEMENT A. Approved for Public Release; Distribution is Unlimited. The authors would also like to thank Stephen Michielsen for the helpful discussions.

Conflict of Interest

All authors declare no conflicts of interest in this paper.