Systematic approach for wettability prediction using molecular dynamics simulations †

We present a fast and eﬃcient approach to predict the wettability and spreading of liquids on polymeric substrates. First, a molecular dynamics parameterization is proposed for the calculation of the solubility parameter for 74 compounds including surfactants typically used in inkjet printing. Then, we introduce a molecular geometrical factor to relate the solubility parameter to the surface tension, obtaining estimates in remarkable agreement with experiments. By using a modified Young–Fowkes equation, the contact angles of liquids on various polymeric substrates are determined and their dependence on the hydrogen bonding, dispersion and polar contribution of the solubility parameter are investigated. We find that wetting properties are obtained with a good accuracy when taking into account the hydrogen-bonding and polar interactions in the geometric sum of the solubility parameter. Based on these findings, a 3D wetting space is proposed to evaluate liquids wettability and judge their suitability for specific substrates. This will enable easy formulation of liquids with wettability tailored for a particular surface and application.


Introduction
Wettability and spreading of liquid drops on solid substrates is a fundamental phenomenon in nature, important for diverse applications, including inkjet printing, 1,2 oil recovery, 3 agriculture 4 and textiles. 5,68][9] As a consequence, optimizing inkjet processes requires formulating ink mixtures that have a specific predefined range of wetting, which is a time-consuming and prohibitively expensive process.Therefore, a rapid and economical approach for the prediction of the wetting properties of liquids is necessary, especially when the screening of thousands of potential mixture components is required.
The wettability of a liquid droplet on a solid surface is determined by the molecular interactions of the materials involved.
Strong attractions between the molecules of the liquid and the solid favors the spreading of the liquid on the surface.Conversely, weak attractions cause the liquid drop to bead off.Under static conditions, wettability is commonly described in terms of the contact angle y (see Fig. 1) by Young's equation, which has been proved recently to be valid down to the nanometer scale, 10,11 cos y ¼ g sv À g sl g lv ; where g sv , g sl and g lv are, respectively, the solid-vapor, solid-liquid, and liquid-vapor interfacial energies.Tuning the wettability requires the balance between these energies to be manipulated.However, difficulties in obtaining reproducible contact angles have been reported by many investigators, and both g sv and g sl cannot be measured reliably. 12,13As a consequence, several semi-empirical approximations have been proposed to predict the contact angle and the surface energy of solids and liquids.These include the critical surface tension approach by Zisman and co-workers, 14,15 harmonic-mean approximation by Wu, 16,17 the acid-base approach by van Oss 18 and the geometric-mean approximation proposed by Girifalco and Good, 19 and by Fowkes. 20Fowkes 20 and later, Owen and Wendt 21 suggested that the surface free energy of a material is the sum of different intermolecular interactions, namely polar, dispersive and hydrogen-bonding interactions.On the other hand, Neumann 22,23 proposed a model based on macroscopic thermodynamics that disregards the concept of surface free energy components.
][31][32] Sgouros et al. 29 developed a simulation strategy for the prediction of the interfacial free energies and wetting properties of polyethylene-graphite interfaces.Aqra et al. 31 proposed a theoretical model based on classical statistical thermodynamics for calculating surface tension of liquid metals.Despite these efforts, there is still little consensus on the determination of interfacial energies, and the separate effect of the disperse, polar and hydrogen-bonding interactions on wettability is still not well understood.
One important thermodynamic property of a material is the Cohesive Energy Density (CED).This property, coined by Hildebrand and Scott as the ''solubility parameter'', 33 is equal to the amount of energy per unit volume required to separate a molecule from its surrounding environment, toward a distance where it approaches zero potential energy.Because it contains all the intermolecular energies holding the liquid together, it represents an important descriptor for liquid wettability.7][38][39][40] Yet, its experimental determination is still difficult, and the solubility parameter of a large number of liquids and miscible mixtures, commonly used in inkjet printing processes, are still unknown.
Typical modeling approaches for the determination of the solubility parameters of molecules are based on theoretical calculations via group contribution methods proposed by Van Krevelen, 41 Hansen, 42,43 Yamamoto (Y-MB) 44,45 and Hoy, 46 where all contributions of the functional groups (e.g.-CH 3 and -OH) are summed up.][49] Molecular Dynamics (MD) simulations is a promising more robust alternative to compute the molecule's solubility parameter.1][52] MD has proven a convenient tool for the calculation of solubility parameters for pharmaceutical products, 53 solvents, 54 CO 2 -expanded mixtures, 55 non-ionic surfactants, 54 ionic liquids, 56 and polymers, 36,[57][58][59] but its use for wettability-related applications has been largely unexplored.
In this paper, we propose a simple and fast approach to predict the wetting properties of liquids.By performing short molecular dynamics simulations, the Hansen solubility parameters 42,43 of molecules are computed more accurately than in methods purely based on the atomic structures of the molecules.The solubility parameter of liquids, surfactants and polymers (74 compounds) obtained from MD are compared against experimental data, showing good agreement.An alternative form of the equation of Beerbower 60 relating surface tension to the solubility parameter is proposed, by including details of the geometry of the molecule, resulting in a better fit to the experimental data.The expression for the contact angle by Owens and Wendt 21 is extended to include all three Hansen solubility parameters, to improve its predictive value.Akin to the miscibility spheres of liquids in the Hansen's 3D space, the revised contact angle expression is represented by 3D wettability spheres that allows to evaluate the performance of liquids in terms of wetting properties on a given substrate.

Solubility parameter and Hansen 3D space
A molecule which ''evaporates'' escapes from the droplet to the vapour phase when its energy becomes high enough to overcome the cohesive energy of the surrounding molecules.This energy of vaporization can be obtained experimentally from the molar enthalpy of evaporation, which, if the vapor phase in equilibrium with the liquid behaves ideally, can be related to the cohesive energy density, whose square root is known as the solubility parameter, 33,61 where U is the molar internal energy of the liquid phase, H vap is the molar enthalpy of vaporization, R is the gas constant, V m the molar volume of the liquid and T is the temperature.Crowley et al. 62 and later Hansen 42 proposed to split the solubility parameter d into three components based on intermolecular interactions.The first component d d is related to the van der Waals dispersion interactions, the second component, d p , is related to the coulombic interactions (excluding hydrogenbonding interactions), and the last one, d h , represents the hydrogen-bonding interactions.The total solubility parameter is the geometric sum of these components, Following this scheme, Hansen 42 developed what is now the most widely used three-dimensional interaction space, wherein all the solid and liquid substances are localized.The basic premise of Hansen's theory is that substances that are closer to each other in Hansen space generally have higher miscibility, see Fig. 2. For a solid substance to be soluble in a liquid, or two liquids to be miscible, it is necessary that their positions in the Hansen space are close, that is to say that their three solubility parameters are similar.Once a wide range of substances are placed in the Hansen space, it will become a useful tool for liquids formulation with less empirical front-end testing.

All-atom molecular dynamics
Molecular dynamics simulations were conducted using LAMMPS. 63The atomic charges of each molecule were obtained using the charge equilibration (Q eq ) method of Rappe ´and Goddard, 64 except for carbonyls (e.g., ethylene carbonate and acetic acid) where the Electrostatic Potential Fitting (ESP) 54,63 was used since it produces better results for these molecules.
The initial amorphous configurations were generated using a Monte-Carlo algorithm available in the Scienomics MAPS package. 65,66The energies of the amorphous configurations were minimized using the steepest descent method, followed by the conjugate gradient method with a maximum number of iterations equal to 1000.After generating the molecular simulation cell with periodic boundaries, DREIDING 67,68 forcefield was applied to the molecules, where 12-6 Lennard-Jones and Coulomb electric potentials account for the van der Waals and electrostatic interactions, respectively, and a 12-10 Lennard-Jones potential accounts for the hydrogen-bonding interactions.This latter was included through a donor-acceptor distance-angle dependent definition. 69,70][73][74] We set the angle cutoff to 1401 to limit the scope of the Lennard-Jones hydrogen-bonding potential, see Fig. 3.All oxygen and nitrogen atoms were considered to be potential hydrogen bond acceptors and donors.The hydrogen bond Lennard-Jones potential was parametrized with a zerocrossing distance R hb = 2.75 Å and a hydrogen-bond energy D hb = 4 kcal mol À1 , with a cosine angle periodicity equaling 4 as suggested by Mayo et al. 67 A cutoff distance of 10 Å smoothed at 11 Å was used for the 12-6 Lennard-Jones and Coulomb interactions.Long-range electrostatic interactions were handled by the Ewald summation, 75 making the Coulomb cutoff distance effectively infinite.The total number of molecules in the simulation cell varies between 100 for the smallest molecules and 10 for the largest ones.All initial configurations were relaxed to adjust the coordinates of the atoms, where the experimental density of the different compounds was used to set the simulation box volume.Molecular dynamics simulations were launched in the canonical ensemble NVT for 200 ps, with a time step of 1 fs.The temperature was set at room temperature T = 298.15K and controlled by a Nose-Hoover thermostat.For polymeric substrates, we followed the same protocol used by Jawalkar et al., 58 Diaz et al. 76 and Jarray et al. 36 Although a degree of polymerization of 8 with four polymer chains are sufficient to compute solubility parameters, 36,54 all our simulations were performed with at least 10 polymer chains in the simulation box with 8 repetition units each.To ensure that the equilibration and sampling durations are adequate, we created independent initial configurations following the same recipe as well as performed longer simulations of 1 ns.][79] Once molecules with promising properties have been identified, they can be subject to longer simulations for more accurate estimates and/or tested experimentally.
The solubility parameter values, with units of [energy density 0.5 ], where calculated by averaging the intermolecular energy 54 over the last 20 ps of the simulation, where k runs over the dispersive, coulombic and hydrogenbonding energies.The brackets denote a time average, and V denotes the simulation box volume.This journal is © The Royal Society of Chemistry 2020

Surface tension and wettability prediction
Surface tension, i.e. the imbalance of the cohesive forces between molecules at the liquid-vapour interface, is an important parameter for assessing the wettability of liquids.It can be related to the solubility parameter using the empirical equation 60,61 where k is a constant equal to 0.0171.Although this equation gives satisfactory values of the surface tension for several liquids, its accuracy can be improved by discerning the types of interactions.Beerbower 60 provided, on the basis of a computer regression of experimental data, a more accurate correlation that differentiates the three components of the solubility parameter, with k the same constant as in eqn ( 5), and b equal to 0.632.Since these constants were derived empirically, it would be nice to derive their physical origin.Becher 80 attributed the constant k to the geometrical properties of the liquid molecule.The difficulty, as seen at that time by Hildebrand and Scott, 61 lies in the evaluation of the surfaces and volumes of the molecules.
Becher assumed that the molecules are highly convoluted and used the characteristic dimensions of the molecule. 80We propose to extend this argument by including the molecule details.Taking into account the three components of the solubility parameters, we define the solubility parameter-surface tension relation as with V a the volume of the molecule, S a the surface area of the molecule, and a a parameter that depends on the ratio of the number of neighboring molecules on the surface to that in the bulk.Roughly speaking, molecules at the surface of a droplet lose about half of their cohesion energy, 81 hence, a E 0.5 as it will be shown in this work.A detailed derivation of eqn (7), explaining the origin of a and b, is presented in Appendix A.
For the computation of V a and S a , the procedure of Richmond 82 and Richard 83,84 was used.The approach consists of rolling a spherical probe on the molecule atoms while in contact with the van der Waals surface, which is defined by the boundary of the union of the atomic van der Waals spheres.The rolling probe center provides the total accessible area of the molecule, see Fig. 4. The probe radius suggested by Richard 83 and often used by many authors [85][86][87][88] is r a = 1.4 Å, and corresponds to the average radius of a water molecule.However, recent studies [89][90][91] suggested the use of a probe that depends on the atom types of the molecule to get better accuracy of the solvent accessible area.Here, we will use a probe radius that corresponds to the average atomistic radii of the molecule, which are provided by Wang and Hou 90 (1.4 Å for oxygen, 1.72 Å for carbon, 1.54 Å for nitrogen, 2.1 Å for silicon, 1.6 Å for fluorine and 2 Å for bromine).As we will show, a constant probe radius of r a = 1.4 Å also gives good results, indicating that 1.4 Å is an acceptable average radius.
Based on Fowkes' approximation 20,92 of additive contribution to the surface energy and Good-Girifalco's geometric mean approximation 93 of interaction energies, Owens and Wendt 21 used the Young equation to postulate the following equation for estimating either the interfacial energy or the contact angle of liquids on solid surfaces: where g s and g l are the interfacial energies of the substrate and liquid, respectively, and the subscripts d and p denote, respectively, the dispersive and polar components of the interfacial energy.g l is the total liquid surface tension given by the sum of g l d and g l p .This geometrical mean approach uses only two categories of interactions instead of three (i.e., polar and dispersive interactions).By taking into account the hydrogenbonding contribution and equating the last three equations, we obtain the wettability in terms of the three solubility parameter components: Eqn ( 9) derived from the empirical eqn (6) uses the molar volume, while eqn (10) derived from eqn (7) uses the molecular geometrical factor ffiffiffiffiffiffiffiffiffiffiffiffiffi V a =S a p and assumes identical values of a for solids and liquids (i.e., a l = a s ).While keeping in mind that both of these equations ignore the effect of surface roughness and inhomogeneities, they can be used to predict y if the various dispersion, polar and hydrogen bonding components of the solubility parameter are known.5][96] Solubility parameter values are summarized in Table 1.Obtaining the three components of the solubility parameters is difficult experimentally, which is why only the total experimental solubility parameters are listed for comparison with the numerical results.A pictorial of the molecular structure of all the compounds used in this study is provided in Section S2 of the ESI.† MD gives solubility parameter values in agreement with experiments with an overall deviation of 1.1%.The largest deviations of 5-10% are observed for amine-based liquids (e.g., formamide) and carbonyls (e.g., ethyl acetate).As previously pointed out by Belmares et al., 54 this discrepancy is because the Lennard-Jones potential of the hydrogen-bonding interactions in the DREIDING forcefield is not well adjusted for nitrogen and double-bonded oxygen atoms.
Fig. 6(a and b) show the experimental surface tension of typical liquids, polymers and surfactants used in inkjet printing, plotted against the surface tension values obtained using eqn ( 6) and (7), respectively.Nonpolar (e.g., hydrocarbons, halogenated -carbons), partially polar (e.g., acetone and ethyl acetate) and hydrogen-bonding compounds (e.g., formamide, glycerol and water) liquids are shown with different symbols.Calculated and experimental surface tension values of these liquids are listed in Table 1.For some compounds showing dispersed experimental surface tensions in the literature (deviation ranges B1 mN m À1 for liquids and B2 mN m À1 for polymers), the average value was taken when plotting Fig. 6.Fitting a and b of eqn (7) to the experimental surface tension values gives b = 0.6225, which is almost equal to the value proposed by Beerbower 60 (b = 0.632), and a value of a = 0.504, showing that molecules on the surface lose about half their cohesive energy.We point out that the value of a depends on the interplay between the number of neighboring molecules on the surface of the liquid and the anisotropic interactions bearing on the surface molecules (see Appendix A).Both equations give surface tension values in good agreement with experiments, but eqn (7) shows an overall better agreement, especially when predicting the interfacial energy of polymeric substrates: In fact, eqn ( 6) and ( 7) agree with the experimental values with 4.3% and 2.4% deviation in comparison to the experimental values, respectively.This endorses the argument of the molecular geometrical origin of the constant k in Hildebrand-Scott 61 and Beerbower 60 equations, previously derived from experimental extrapolation.The agreement also confirms that better predictive accuracy is obtained when taking into account the type of interaction, namely polar, disperse and hydrogen-bonding.In particular, the surface tension of hydrocarbons is well reproduced and fitting eqn (7) to only hydrocarbons gives b = 0.96 and a = 0.52 (see Appendix B).Alkane molecules are dominantly dispersive, and have weak orientational ordering due to their weak polar and hydrogenbonding interactions, which explains why b has a low effect on their surface tension.][105] From this standpoint, we can infer that the constant b in the Beerbower 60   This journal is © The Royal Society of Chemistry 2020 Table 1 Experimental and computed solubility parameters in (J cm À3 ) 0.5 obtained using MD simulations, volume (Å 3 ) and surface accessible area (Å 2 ), experimental and calculated surface tensions (mN m À1 ) at 298.15 K, obtained using eqn ( 6) and ( 7 hydrogen-bonding interactions lost by the similarly oriented molecules at the vapour-liquid interface.When using a constant probe radius r a = 1.4 Å for the calculation of the solvent accessible area in eqn (7), good results were also obtained as shown in the ESI † (Section S3), with a 3% deviation in comparison to the experimental values.Solubility parameter and surface tension values obtained using the Yamamoto-Molecular Break (Y-MB) method 44,45 are also provided in the ESI † (Section S4).These results show that the Y-MB method works equally well for small molecules such as glycerol and cyclobutanone, but becomes less accurate for polymers and large molecules such as tergitol.Overall, better results are obtained when using the solubility parameters components computed from MD.

equation takes into account the polar and
Next, we explore the wetting properties of liquids on different smooth polymeric substrates namely PMMA: polymethyl methacrylate, PTFE: polytetrafluoroethylene, PP: polypropylene, PC: polycarbonate, PET: polyethylene terephthalate, POM: polyoxymethylene, PS: polystyrene, and PVA: polyvinyl alcohol.In Fig. 7, the calculated contact angles obtained using eqn (10) are plotted against experimental data from literature, obtained using the sessile drop method and the one-liquid contact angle method. 100,107All contact angle data are tabulated in the ESI † along with a version of Fig. 7 where liquids are labeled with their name (Sections S5 and S6, respectively, ESI †).The polymeric substrates are smooth with a surface roughness factors below 1.05, i.e. the ratio of the actual surface area of the substrate to its horizontal projection, thereby minimizing the estimated roughness-induced enhancement of the experimental contact angle value to less than 2 degrees. 108The calculated and reported experimental contact angles are equilibrium values, which usually lie in between the receding and advancing angles of moving droplets. 109,110Our equation, with b = 0.6225, gives results in reasonable agreement with experiments, which favors Good's 93 geometric mean in eqn (8).
The Teflon (PTFE) substrate, a well-known low-surface-energy materials that is frequently used in fusing components in the printing industry, 111,112 shows the best agreement.Similarly, the computed contact angle values of polar and H-bonding polymers, such as PC and PET, are close to the experimental values.However, for PS and PP, there is an overestimation of the contact angle when they are in contact with highly polar liquids, but the equation correctly predicts the expected trend.This may have several causes: for instance, on the experimental side, beside the difficulty of the acquisition of precise contact angle values, the effects of surface heterogeneity can easily overshadow the influence of interfacial interactions. 113,114On the theoretical side, as stated by Wu 102 the geometric mean in Abbreviation: PC: polycarbonate, PET: polyethylene terephthalate, DPGME: di(propylene glycol)methyl ether, PMMA: polymethyl methacrylate, POM: polyoxymethylene, PP: polypropylene, PTFE: polytetrafluoroethylene, PS: polystyrene, and PVA: polyvinyl alcohol.For the computed solubility parameters, the standard deviation is in the range of B0.2 (J cm À3 ) 0.5 for small molecules and B0.25 (J cm À3 ) 0.  This journal is © The Royal Society of Chemistry 2020 eqn (8), sometimes does not perform well for interactions between low energy materials (e.g., hydrocarbons) and highly polar liquids (e.g., water and glycerol), especially if they involve molecular rearrangement at the solid-liquid interface resulting from weak H-bonds formed between the hydrogen atoms of the substrate and the hydrogen acceptors of the liquids.In this particular type of overwhelmingly dispersive hydrocarbon substrates, we recommend modeling their interactions with polar liquids on the basis of the total solubility parameter rather than with their individual separate interaction components, which curiously gave slightly better results for PS and PP (and for these two only), as shown in the ESI † (Section S7).Alternatively, considering the orientation of the molecules at the solid-liquid interface and taking into account the dependence of a and b on the type of the interface in eqn (10) (i.e., taking a s a a l rather than assuming a constant a) may provide more accurate predictions.
To enable numerical vs. experimental comparison of the spreading behaviour of the liquids, we will consider next that liquids completely spread on the substrate if the calculated contact angle value is lower than 151.The wetting behaviour is summarized in Fig. 8. Bottom-right and upper-left oriented triangles denote experiments and numerical calculations, respectively.Red and blue triangles denote dewetting and spreading liquids, respectively.To interpret the significance of the contact angle results, components are sorted according to the magnitude of their hydrogen-bonding interactions.Wetting behaviour computed numerically agrees with experiments in all cases.According to Fig. 8, strongly hydrogen-bonding liquids dewet or partially wet all the polymeric substrates.PTFE is neither a hydrogen bonding nor a polar substrate, which explains why no liquids spreads on it.Even for the strongly hydrogen-bonding liquids, polar and dispersion interactions still play a role.For instance, the latter explains why bromonaphthalene partially wets PMMA but spreads on the less polar PC surface, in agreement with experiments (see Table S3 in the ESI †).
Similarly to the Hansen 3D space approach, the three parameters d d , d p and d h , multiplied by the square root of the ratio of the solvent volume to the accessible solvent area ffiffiffiffiffiffiffiffiffiffiffiffiffi V a =S a p , can be treated as coordinates of a point in threedimensional space.Then, by solving eqn (10) for a specific substrate and a predefined contact angle, a wettability sphere is obtained as shown in Fig. 9.Note that the coefficient b of eqn ( 10) is taken into account in the coordinates of the polar and hydrogen-bonding components to generate a sphere instead of a spheroid.Here again, in eqn (10), we consider b to be a constant equal to 0.6255.Liquids that lie inside the resulting sphere for a certain predefined low contact angle, will spread on the corresponding substrate.
In Fig. 9, we plot two wetting blue and red spheres of PTFE (Teflon) corresponding to contact angles of 151 and 901, respectively, and we place the coordinates of the 66 liquids investigated in this study, whose properties are tabulated in Table 1.The blue and red dots in Fig. 9 present, respectively, liquids which are inside the blue and the red wetting spheres.The blue dots, inside the blue sphere, mark liquids that spread; the red dots, in between the red and blue spheres, indicate liquids that  form less well spread droplets up to hemi-spherical droplets, while the black dots represent liquids that dewet the substrate (i.e., form a spherical droplets on the substrate).We notice that as we go further from the origin of the 3D wetting space, liquids spreading decreases because the polar, dispersive and hydrogen-bonding contributions to the cohesive force between the molecules of the liquid increases.In the specific case of PTFE, almost all liquids lie outside the blue wetting sphere, indicating dewetting or partial wetting.Only a few hydrocarbon liquids (e.g., toluene) are near the edge of the blue sphere and will likely spread.The other liquids, in particular water and glycerol, exhibit a high polar contribution hindering the wetting of PTFE.We also provide the wetting spheres of all the other polymers in the ESI † (Sections S8).From Fig. 9(c), it is interesting to see that for many liquids the hydrogen-bonding component increases almost linearly with the polar components of the solubility parameter.
Knowledge of the wetting sphere allows easy wettability determination and fast screening of adequate solvents and surfactants for many applications including inkjet printing.Here for instance, the wetting sphere of PTFE needs to be enlarged to improve the spreading, for example, by means of surface treatment, or in the case of water and glycerol, by selecting a surfactant whose polar contribution is low.Furthermore, this approach can also be applied to a combination of miscible liquids, since MD allows the computation of the solubility parameter of miscible mixtures.For rough surfaces, eqn (10) can be extended by using for example the Wenzel model, 108 simply by multiplying its right side by the roughness factor.

Conclusions
We introduced a straightforward, low-cost, molecular dynamicsbased approach for characterizing the wetting behavior of liquids, speeding up the screening of adequate molecules that meet predefined wettability and spreading requirements encountered in various industrial applications.We proposed an MD framework that predicts well the solubility parameter components of popular liquids and polymeric substrates.Then, we showed that the empirically derived Hildebrand-Scott constant in Beerbower's correlation relating the surface tension to the solubility parameter, 60 is based on geometrical properties of the liquid molecule.By separating the solubility parameter components in Owens and Wendt 21 equation and including the hydrogen-bonding interactions, we proposed a new 3D wettability space that is inspired by the Hansen solubility parameter space, 42 to generate the wetting spheres of substrates.This allows to categorize the wetting behaviour of liquids or mixtures and select or design suitable substrates for specific applications.
The proposed approach is fast and yields Hansen solubility parameters, surface tensions and contact angles in good agreements with experimental data.One may therefore expect the method to give reasonable predictions for compounds that were not part of the current study.Compared to group contribution methods such as Y-MB, 44 the method is not as fast but it is more accurate, accounts for conformation changes in molecular structure, and can be applied to polymers and large molecules.In comparison to studies extracting Hansen parameters, surface tensions or contact angles by simulating long polymer chains or explicit droplets, 29,[115][116][117][118] the current method is much faster.There are a few limitations that need to be addressed in future studies.For instance, the DREIDING 67 hydrogen-bonding potential is still not well parametrized for nitrogen and bromine atoms, and can be adjusted for higher accuracy.Additionally, Good's 93 geometrical mean approximation implemented in our equation sometimes overestimates the contact angles between dominantly dispersive substrates and highly polar liquids.Consequently, it would be interesting to further universalize our equation to accurately predict liquid wettability while taking into account these specific interactions.
We want to determine the physical origin of the constants a and b in eqn (7).The surface energy of a pure liquid is the energy needed to increase the interfacial area by one unit.It can therefore be approximated as the energy required to bring one individual molecule from the bulk to the surface of the liquid, 119  where e is the average pairwise intermolecular energy, S is the surface area of the molecule exposed to the vapour phase, and Z is the coordination number.The denominator ''2'' is due to the double-counting of the pairwise interactions, and the superscripts s and b refer to the surface and bulk of the liquid, respectively.Similar equations were proposed by several authors. 31,120,121However, these authors assumed that the molecules on the surface interact the same way as in the bulk (i.e., e s = e b ).This should not be the case, particularly for polar and hydrogen-bonding molecules as their orientation may change on the vapour-liquid interface, subsequently affecting the intermolecular interactions undergone by the molecule on the surface.Assuming that the intermolecular interactions are equally distributed over the surface of the molecule, the surface area exposed to the vapour phase can be written in terms of the total accessible area of the molecule as S = S a (Z b À Z s )/Z b .Introducing the ratio b k = e s k /e b k , where k runs over the van der Waals dispersion, the coulombic, and the hydrogen-bonding energies, and using we get The ratio b k takes into account the excess binding energy.This latter equation translates into eqn (7), when taking Here, a depends on both the number of neighboring molecules as well as on the intermolecular energy on the bulk and on the surface.Note that this derivation of the surface tension as a function of the solubility parameter is energy-based and disregard the effect of the entropy, which is taken into account in the fitting of a and b in eqn (7) against experimental data.

Fig. 1
Fig. 1 Schematic of the wetting of a droplet on a substrate (right).High wetting is characterized by a low contact angle.Wetting is determined by the interactions at the atomistic scale (left).

Fig. 3
Fig. 3 Geometrical MD parameters for describing hydrogen-bonding between two formamide molecules: f AHD is the hydrogen bond angle.Here, nitrogen (in blue) is the donor atom and oxygen (in red) is the acceptor atom.Fig. 2 3D representation of the Hansen solubility parameters.Water and glycerol are miscible because their positions are close in this space.Toluene, on the other hand, has a lower affinity for water because they are far from each other.

Fig. 5
Fig. 5 MD-calculated versus experimental solubility parameters.Experimental values of the solubility parameter are gathered from ref. 94-97.
Abbreviation: PC: polycarbonate, PET: polyethylene terephthalate, DPGME: di(propylene glycol)methyl ether, PMMA: polymethyl methacrylate, POM: polyoxymethylene, PP: polypropylene, PTFE: polytetrafluoroethylene, PS: polystyrene, and PVA: polyvinyl alcohol.For the computed solubility parameters, the standard deviation is in the range of B0.2 (J cm À3 ) 0.5 for small molecules and B0.25 (J cm À3 ) 0.5 for polymers, see Section S1 of the ESI.The corresponding deviation in the surface tension is in the range of B0.6 mN m À1 for small molecules and B0.8 mN m À1 for polymers.a Calculated from experimental enthalpy of vaporization, DH vap from ref. 95 and 98, using eqn (2).b Experimental surface tensions at 293.15 K. c Dispersed experimental values in the literature.d Experimental solubility parameters are obtained from ref. 94-97.e Experimental surface tension values are obtained from ref. 94 and 96-102.f Eqn (7) with a probe radius that depends on the atom types of the molecules.

Fig. 9
Fig.9(a) Wettability spheres for PTFE, (b) projection on the hydrogenbonding-dispersive plan, and (c) projection on the-hydrogen-bondingpolar plan.All data are in (J cm À2 ) 1/2 .Liquid dots inside the blue sphere spread on PTFE, while those outside but still inside the red sphere partially spread on PTFE.Dewetting happens for liquids outside the red sphere.

Fig. 10
Fig. 10 Schematic representation of a molecule on the surface and in the bulk of the liquid.Polar molecules do not interact the same way on the surface as in the bulk because they exhibit a preferred orientation on the surface. )